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ABSTRACT 


The buckling load of a blade stiffened laminated composite plate having midplane 
symmetry 1s maximized for a given total weight. The thicknesses of the layers and the 
width and height of the stiffener are taken as the design variables. Buckling analysis is 
carried out using a finite element method. The optimization problem is solved using 
commercially available optimization packages. 

Due to the highly nonlinear nature of the optimality equations, several local opti- 
mum solutions are found. To examine the relationship between the number of local 
optimums and their relative magnitudes, various combinations of fiber orientation for 


the laminate layers and the blade stiffener are investigated. 
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THESIS DISCLAIMER 


The reader is cautioned that computer programs developed in this research may not 
have been exercised for all cases of interest. While every effort has been made, within 
the ume available, to ensure that the programs are free of computational and logic er- 


rors, they cannot be considered validated. Any application of these programs without 
additional verification is at the risk of the user. 
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I. INTRODUCTION 


Design optimization is not a new subject. Man has always strived to build the best. 
whether it be ancient man making a hunting bow or today’s engineer designing a light- 
weight, high-strength wing for a high-performance aircraft. The goal for both 1s exactly 
the same, to make the best possible design subject to a set of design requirements. The 
difference between the two is not in the goal but in the technique. Ancient man used 
trial and error to improve on his previous designs and, over a period of time, he was able 
to produce an efficient weapon. Today, high-speed digital computers are used to analyze 
and improve designs by solving a large set of equations that have been formulated to 
mathematically model the design problem. 

The design optimization problem of current interest is the maximization of the 
buckling load for a blade stiffened composite plate. This problem can be viewed from 
two different perspectives. On one hand, this is a design optimization problem which 
optimizes the shape of the plate and stiffener cross section to give the maximum load. 
On the other hand, this is a design optimization problem which optimizes the thickness 
and/or the ply angle of the composite laminae to give the maximum load. Previous 
works on the optimization of the plate and stiffener cross sectional shape have focused 
on the use of solid elastic materials [Refs. 1-4]. In the area of design optimization of 
composite plates, previous works have focused on plates without stiffeners. Optimiza- 
tion of the composite plate with respect to the lamina fiber orientations has been studied 
in [Refs. 5-17]. Of greater interest to this study are the previous works on the design 
optimization of composite plates where the laminae thickness and not the ply onen- 
tations are taken as the design vanables [Refs. 18-21}. The current study combines the 
design optimization problems of cross section optimization and composite laminae 
thickness optimization into a single design optimization problem. This is done by opti- 
mizing a blade stiffened composite plate for the maximum buckling load by using the 
laminae thicknesses and the stiffener height and width as the design variables. 

Composites are used as the structural material in this study because of their high 
strength-to-weight and stiffness-to-weight ratios. These properties are very important 
in many of today’s engineered structures, especially in the aerospace industry where high 


strength is required but a severe penalty is incurred for added weight. Because of this, 


composites make an excellent choice for a structural material. However, the use of 
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composites greatly complicates the design analysis problem arising from the use of lam- 
inated materials with non-homogeneous properties. Now the analysis process must 
consider the orthotropic properties and the ply orientation of each individual lamina, the 
specific stacking order of the laminae, and the thickness of each individual lamina. One 
serious problem that arises is the coupling effect between extension and bending. Be- 
cause a symmetric composite plate is chosen, the coupling between extension and 
bending is eliminated. This lack of coupling allows the use of linear buckling analvsis 
up to the buckling load. 

The formulation of a suitable design analysis technique is the first step in the design 
opumization. Since the optimum design can never be better than the design analvsis 
that it is based on, the formulation is described in some detail. The buckling analysis 
of the stiffened plate is performed in two steps. First, a finite element method is used 
to determine the overall buckling load of the stiffened plate. This is done by treating the 
stiffener as a beam in forming the global stiffness matrix for the plate. Using FEM 
formulation, the first three buckling loads are solved numerically. The second part of 
the buckling analysis finds the local buckling load of the stiffener. By treating the 
stiffener as a plate with three sides simply supported and the fourth side as free, the local 
buckling problem can be solved analytically using a Levy's solution. 

Once the design analysis technique has been formulated, the design variables, taken 
to be the laminae thickness and the stiffener height and width, are optimized using an 
optimization program. This study uses three different optimization programs which are 
all commercially available. They are: the IMSL Subroutine DNCONG; the Design 
Optimization Tools (DOT) program using Modified Feasible Direction; and the Design 
Optimization Tools (DOT) program using Sequential Linear Programming. 

The specific problem in. :stigated in this study is a symmetric blade stiffened com- 
posite plate that has two laminae and a stiffener on either side of the centerline. The 
fiber onentation of the laminae and the stiffener are allowed to be either 0° or 90°. Four 
different stacking orders of these fiber orientations are investigated. For each case, de- 
sign optima are investigated. These optima are then compared to determine if there is 
a best stacking configuration. The study is completed for two different total volumes 
of material. Finally, the relative efficiencies of the three optimization routines are in- 
vestigated to determine qualitatively which is the best suited for this specific problem. 
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II. DESIGN APPROACH 


A. GEOMETRY 

The geometry used in this study is a blade stiffened composite plate which is sym- 
metric to the plate centerline as shown in Figure 1 on page 4. The plate is simply sup- 
ported on all four sides and subject to a compressive inplane load in the stiffener 
direction. The blade stiffener is centered on the plate and spans the plate length, L,. 

Each lamina has a thickness of T, and a ply orientation of angle 6,. The stiffener 
is Gefined by a width of B and a height of H. The ply orientation of the stiffener can 
be either 0° or 90° . The independent engineering properties of each lamina and the 
stiffener are given by E),, Ey». Vi2, Giz, and E,,,, Ez, Yi2,, Gia, » respectively. 


B. STRESS ANALYSIS 

The analysis of a composite plate structure is based on classical lamination theory. 
From this theory, the stiffness of a composite structure can be determined from the 
material properties of the individual laminae [Ref. 22: p. 147]. -The plate in this study is 
composed of individual layers which are made of orthotropic material. Also, the plate 
is thin with respect to its span length so it is assumed to be under plane stress conditions. 
From this, classical lamination theory can be applied to determine the composite’s 
stiffness. 

First, the stress strain relation for a orthotropic material under plane stress must be 
defined for the individual laminae. This relationship is given in Eq. (1) below. 


oy Qi Qi2 Me e} 
8 P=] Qi2 Qn Qe 7) (1) 
O12 Qi6 Qr Ve Y12 


where the reduced stiffness matrix elements, Q,, are defined in terms of the elastic 


y? 
modulus in the 1-1 principal direction, £,,, the elastic modulus in the 2-2 principal di- 
rection, £,, the shear modulus, G,,, and poisson’s ratio relating transverse strain in the 
1 direction when stressed in the 2 direction, v,, , and poisson’s ratio relating transverse 


strain in the 2 direction when stressed in the 1 direction, v,, . 
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Figure 1. Problem Geometry: top) general -view bottom) cross sectional 


view 
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For each lamina, the principal axes lie along the direction of the fiber orientation. 
Since the fiber orientation of each lamina is allowed to vary by an angle 6, Eq. (1) must 
be transformed to determine the stress strain relationship in the lamina geometric refer- 


ence frame. This relationship is given in Eq. (2). 


Ox oe Pi Or ey 
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where the elements of the reduced transformed stiffness matrix, 0 
of Q, and @ as 


, » are defined in terms 


Q1, = Q), cos + 2°Q\. + 2Qee) sin’@ cos’?@ + Q,, sin*d 

Orn = (Qi: + Qr2 — 406) sin’@ cos’@ + Q,( sin"@ + cos*@) 

Q4.= Ai, sin'@ + 2(Q12 + 2066) sin’@ cos’@ + Q,, cos"@ 

Ore = (Qi; — Diz — 26) sin 9 cos’ + (Qy2 — Qe. + 206) sin’@ cos 6 
Qo = (Qu — Qi2 — 2Q%66) sin’@ cos @ + (Qi. — Qe. + 2Q%6) sin 8 cos°@ 
Dos = (Qiy + Qrr — 2012 — 26} sin’6 cos + Qgg(sin’@ + cos‘) 


Now Eq. (2) can be thought of as the stress strain relationship for the k® layer of a 


multilayered laminate. Thus it is now written as Eq. (3). 











fo}, = (Q]alehs (3) 


The strain at any point through the depth of the laminate can be expressed in terms of 
the midplane strain, e’, the distance z from the centerline, and the midplane curvature, 
«x. See Figure 2 on page 7 for the definition of z . 


The strain is now expressed as Eq. (4). 
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and 
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The resultant forces and moments acting on a laminate are obtained by integration 








of the stresses in each layer through the laminate thickness. 
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Figure 2. Half of Symmetric Laminate Plate 
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Substituting Eq. (4) into Eqs. (Sa) and (5b) gives, 
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Dy = &D Gye (et - 2h) 
k=l 
The A, terms are the extensional stiffness, the B, terms are the coupling stiffness, and 
the D, terms are the bending stiffness. Since the composite piaie modeled in the study 
is symmetric about the middle surface, the terms B, all drop out [Ref. 22: p. 164). 
Therefore, in this problem there is no coupling between extension and bending. 


C. BUCKLING ANALYSIS 

Buckling of the stiffened plate can occur in two ways. These are: the plate and 
stiffener can buckle as a unit into one of the global buckling modes; and the blade 
stiffener can experience local buckling. Because of this, the buckling analysis is per- 
formed in two steps. For the first case, the blade stiffener is treated as a beam and in- 
cluded in the formulation of the global stiffness matrix in the Finite Element 
Formulation. For the second case, the blade stiffener is treated as a plate subject to a 
inplane compressive load that has three sides simply supported and the fourth side free. 
{Ref. 23: pp. 2-3] 

1. Overall Plate Buckling Analysis 

From [Ref. 23: p. 4], the overall buckling equation in matrix form is 


[K]HU) + CKI(U) — mlKel(U} — polKole(U) = 0 


where [K] and [K], are the plate and beam stiffness matrices, respectively, [K;] and 
[K.], are the plate and the beam geometric stiffness matrices, respectively, and {U} is the 


Nee 


buckling shape. The amount of load carried by the beam and plate are p, and n,, re- 
spectively. The load share carried by the beam and by the plate is proportional to the 


relative stiffness of each. These relative loads are given by 
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where s, and s, are the nondimensional cross-sectional area of the plate and beam, re- 
spectively, and e, and e, are the nondimensional elastic moduli of the laminated plate 
and the beam, respectively. The nondimensional elastic modulus of the plate must be 
determined from the properties of the individual laminae. This value is calculated from 


the extensional stiffness matrix and is given by (Ref. 23: p. 4] as 


ILa]] 


(427466 = 36) tr 


Now the overall buckling equation can be written in the form, 
CK]AU} - p[Kg]f{U} = 0 


where [K], is the total stiffness of the plate and stiffener combination and [Kg], is the 
total geometric stiffness of the plate and stiffener combination. This equation can easily 
be recognized as an eigenvalue problem with p representing the eigenvalues and {U} re- 
presenting the eigenvectors. This problem is solved using DNLASO, one of the sub- 
routines from the package LASO2 [Ref. 24], which computes a few eigenvalues and the 
associated eigenvectors of a large (sparse) symmetric matrix using the Lanczos algorithm 
[Ref. 25]. The programming code used to calculate the first three buckling loads of the 
suffened plate is Subroutine EIGEN contained in Appendix A. 
a. FEM Formulation 

In order to determine the stiffness matrices for the buckling equation, a fi- 
nite element formulation is required. For this study, a 16 degree of freedom plate ele- 
ment is chosen for the plate and a 4 degree of freedom beam element is chosen for the 
stiffener. The plate is divided into four elements (2 X 2) and the stiffener is divided into 


two elements. 





Fides 





From (Ref. 26: p. [16] the beam element stiffness matrix, [4],, is defined as 


12-6, -12 -6l, 
ei, | -6, 42 6f 22 


k= —> 
ne po | -12 6, 12 64 
-6, 22 6 40 


where the element length, /,, is one half the nondimensionalized span length, /,, and e, 1s 
the nondimensionalized elastic modulus of the stiffener. The nondimensionalized mo- 


ment of inertia of the stiffener, i, is defined as 
: 2 3 3 
lb = + oft + ty - 1] 


where ¢; is one half 1¢ total plate thickness for the symmetric plate. From [Ref. 26: p. 


388] the beam element geometni stiffness matrix, [k,],, is defined as 


Free 


a3). 42> 3b. =12 
—36 31, 36 3, 
~3, -12 3 48 


Po 
[ke], = 301 
e 


For the plate element, the geometric stiffness matrix elements, ks, , are 
given by (Ref. 27: p. 432] as 


1th’ af of 
kg, = nf I (SE )( lace (7) 


where f are the 16 shape functions for the element. The shape functions, f, are given 


below, 


fi = (1 — 36? + 26)(1 — 39? + 2n°) 


fy = (1 — 367 + 26°)(n — 2n? + 11°) 
fy = -(€ -— 26 + YL - 3? + 2’) 
fa = (E - 267 + EY(n — dn? + 7’) 


10 


fy = (36? — 28)(1 - 39? + 27°) 
f = (36° ~ 26°)(n — 21? + 1’) 
fr = (2 = &)\(1 — 3n? + 2n”) 

f= —(F - &Y(n ~ 27 + 2’) 


fy = (38? ~ 28°)(3y? — 27°) 


il 


a _ — (36 = 22?)(n? -_ n°) 


il 
— 
wv 

Nv 


wie ~ &)(3n° — 2y’) 
fro = (2 - &)( - 9) 


fig = (1 — 32? + 26°)(3y? — 27) 


fia = —(1 — 36? + 28)(n? - 9°) 
va = ~(é = 22? + &*) (30° _ 2n’) 
fie = ~(E — 22 + &)(w - 9°) 


The integral of Eq. (7) is solved numerically by using gauss quadrature. 
For the plate element, the stiffness matrix elements, k,, are given by (Ref. 
27: pp. 415-418] as 


il 
kj = | | { DC, + DC, + DiC + DC, + DyCs a DeeCo}dédn 
0 “0 


where 





reve 








a2 
Cf, 
Cue ( Sean i Ee 


Again, here the integration is done numerically by gauss quadrature. All of the ele- 
mental suffess formulation is done in Subroutine ELESTF in Appendix A. 

Now that all the elemental stiffness matrices (both stiffness and geometric) 
have been determined, the individual elemental matrices can be assembled into the global 
stiffness matrices. The geometric elemental stiffness matrices [kK,] and [k,], are com- 
bined to form the total global geometric stiffiess matrix [K,];. The elemental stiffness 
matrices (k] and [k], are combined to form the total global stiffness matrix (K],. This 
global stiffness matrix assembling is done in Subroutine ASSEMB in Appendix A. 

2. Stiffener Local Buckling Analysis ; 

To calculate the local buckling load of the blade stiffener, it is treated as a sep- 
arate plate buckling problem. Now the stiffener alone is treated as a plate with three 
sides (x = 0, x = 1, z = 0) simply supported and the fourth side (z = A) free. It is 
subject to axial compressive load p,. Because of the geomeiry and the fact that only the 


first buckling Ic id is of interest, this problem can be solved analytically. 








The general partial differential equation for the plate buckling problem is given 
from [Ref. 22: p. 260] as 
aw aw ow aw 


en te AD 5 + 2D 5) a5: a 5 
Cx Ez Cz Cx 








where now the D, refer to the properties of the blade stuffener and p refers to the local 
buckling load p,. By applying the boundary conditions and assuming a Levy's solution 


of the form 


the buckling problem can be solved directly [Ref. 28: p. 208]. The calculation of p, 1s 


carned out in Subroutine LEVYP contained in Appendix A. 


frarey 











Il. OPTIMIZATION PROBLEM 


A. PROBLEM STATEMENT 

The optimization problem is to maximize the buckling load of the blade suffened 
piate for a given total material volume (which is related to the total weight). The design 
variables are set as the nondimensionalized thickness of the individual lamina, ¢, and the 
nondimensionalized width and height of the stiffener. Here the nondimensionalized 
width, 6, and the nondimensionalized height, A, of the stiffener are denoted as ¢,., and 
t..,, respectively. The nondimensional design variables, 1, are subject to side constraints 
of 


F SG SU may for = 1,2,...,(4 + 2) 


where #,,,, and £,,,, are upper and lower bounds on the design variable dimensions, re- 
spectively, and n is half the number of layers for the symmetric plate. 
The optimization problem for maximizing the buckling load is wnitten as 


i B 
i 
subject to Scr, + 210:5,,-O=0 
1'000p, >B 
0.999p, = B 
0.998p, > B 
p,=B 


and lian SESH for i = 1,2,...,.(2 + 2) 


where p,, p, and p, are the first three overall buckling loads of the stiffened plate, p, is the 
local buckling load of the stiffener, ¢, are the design variables, © is the 
nondimensionalized total plate volume, and c is the relative width of the plate L,/L,. £ 
is a parameter introduced to raise all the buckling loads during the optimization process, 
which finally converges to the buckling load (lowest eigenvalue) of the optimized design. 


The first constraint represents the total volume constraint. Simply put, the total 


volume of the design cannot exceed the maximum available resource. The next four 








ergt 


constraints are minimum limits on the buckling loads. Here the buckling load of the 
stiffened plate must always be at least as large as the smallest value of p,, p, , p; OF P,. 
The coefficients of 1.000, 0.999, and 0.998 in the second, third, and fourth constraints. 
respectively, are necessary to allow the calculation of the eigenvalues, p, , p,, and p,. 
when the buckling mode is bimodal or trimodal [Ref. 29: p. 355]. Finally, the side con- 


straints place minimum and maximum limits on the size of each design variable. 


B. OPTIMIZATION METHODS 
The optimization problem above can be solved by a number of different methods. 

For this study three different commercially available programs are used to solve the op- 
timization problem. Each of these methods uses a different optimization technique 
thereby providing three independent solutions of the optimization problem. The three 
methods used are: the IMSL library Subroutine DNCONG; Design Optimization Tools 
(DOT) using Modified Feasible Direction; and Design Optimization Tools (DOT) using 
Sequential Linear Programming. Each of these methods are variations on the general 
optimization routine. The general steps in the formulation of an optimization routine 
are listed below. 

1. Start with initial guess X° for the 0% iteration, g = 0. 

2. Update the iteration number, g = q + I. 


3. Evaluate the objective function, F(X), and the constraint functions, g(X), at the 
current value of X. 


4. Identify the critical constraints, J. 


5. Calculate the objective function gradient, VF(.X), and the critical constraint func- 
tion gradients, Vg(.X), for je J. 


6. Determine the search direction, S*. 
7. Perform a one-dimensional search to find the step length, a’. 


8. Update the solution by adding the previous iterate solution to the product of the 
step length and search direction, X*' = X™' + a°S*. 


9. Check for convergence. YES: done; No: go to step 2. 


For each of the optimization routines used in this study, a description of the variations 
from the standard optimization routine are given below. 
1. IMSL Subroutine DNCONG 
DNCONG is based on Subroutine NLPQL which is an optimization program 
developed by Schittkowski. This programming algonthm uses a successive quadratic 


programming method to solve the general nonlinear programming problem. In this 


i 


agit 








method, the search direction subproblem is formulated and solved by using a quadratic 


approximation of the Lagrangian function and by linearizing the constraints. The 


Lagrangian function, L(X, 2), is given by 


LIX, 3) = AR — YA9(X) 
j=! 


where / is the vector of Lagrange multipliers and m’ is the number of function con- 
straints plus side constraints, m’ = m+2n. The search direction subproblem is for- 
mulated as 

min 


mae 1 sTpis + vixt)’s 
Ee 2 


subject to Ve(X7)'S + o(X7) = 0, f= lyme 
Ve(X4)'S + g({X7) 20, j=mtl,..m 


and Xan —~ X'S SS XA for k=1,2,...,” 


where 8 is a positive definite approximation of the Hessian matrix of the Lagrangian 
function and X? is the current iterate. The solution of the subproblem, S*, is used as the 
search direction in the line search to find the new point X*'. (Ref. 30] 
2. DOT using MFD 

Design Optimization Tools (DOT) by VMA Engineering is a programming code 
designed to solve a variety of nonlinear constrained or unconstrained optimization 
problems. The version of DOT that uses a Modified Feasible Direction algorithm has 
a search direction, S*, that is a modified form of the steepest descent direction. The ac- 
tual modification to the search direction ts done by using the Fletcher-Reeves conjugate 


direction method. With this method, the search direction is defined as 
Si = —VF(XT™') + CST" 


where 


nlzitsa lk 
vA) P 





reve 


This is a first order search direction which is extremely simple to calculate. The modifi- 


cation of the search direction gives dramatically improved results over the method of 


steepest descent. (Ref. 31: p. E-16] 
3. DOT using SLP 

Design Optimization Tools (DOT) using Sequential Linear Programming 1s a 
method that approximates a nonlinear problem as a linear problem and then optimizes. 
First, a first order Taylor Series approximation of the objective and constraint functions 
are calculated. Then, this approximation is optimized instead of the original nonlinear 
functions. Since the problem is now linear, the value of the objective and constraint 
functions are easily and inexpensively calculated. Also, now the gradients of the objec- 
tive and constraint functions are available directly from the Taylor Series expression. 
The linear approximation problem is optimized using the method of Modified Feasible 
Direction. (Ref: 31: p. E-33}] 


etyt 








IV. RESULTS AND DISCUSSION 


The laminate plate chosen in this study has two lavers with equal span length, L.. 
and span width, L,. The material selected for the laminae and the stiffener is a 


graphite epoxy composite. The material properties of this composite are: 


Ey, = 31.0x 10° psi, Ey, = 34x 10° psi 


G1, = 0.75 x 10° psi, v;7 = 0.28 


Design optimization results are obtained for two different total volume conditions (O 
= 0.04 and 0.015) and four different fiber orientation stacking sequences: 


(0° cam | 90°/ 0°) sym 


OO beam / 0°/ 1) ae 


ras 


(0 aeani | 90°] og es 
(90° seam / 0°/ 90") im 


The minimum and maximum thickness of each laminae are limited to 0.001Z, and 
0.2L,, respectively. The same upper and lower limits are used for the stiffener width and 


height. 


A. OPTIMUM SOLUTIONS 

For each configuration, the known optimums and the values of the design variables 
at these optimums are given. Also listed are the active buckling modes and the number 
of the corresponding figure which presents the same information graphically. All of the 


results are presented in non-dimensional form. 








1. Case 1: © = 0.04 and Ply Orientation (0°... / 90°/ 0°),,. 


For this configuration, two different optimum solutions are identified. They are 
listed in Table 1 below. The same results are presented graphically in Figures 3-6. 


Table 1. RESULTS: @ = 0.04, PLY ANGLE (0°yzex/ 90°/ 0°) sag 


Buckling 
ae a ee ae 


0.00738 0.01020 0.03801 0.06369 | 0.0002616 | 1.2.3.4 
0.01521 0.00100 0.05794 0.06526 | 0.0002822 





puget 





Figure 3. 


Geometry for Case | with p = 0.0002646: For Case 1, total volume O 


= 0.04, with ply angle orientation (0°,,,,./ 90°/ 0°), (to scale) 
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Figure 4. Buckling Shapes for Case 1 with p = 0.0002646: For Case 1}, total 


volume © = 0.04, with ply angle orientation of (0°,.,,,/ 90°/ 0°),,,. Re- 
presented above are the buckling shapes of the simultaneous overall 


plate buckling modes. 
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Figure 5. Geometry for Case 1! with p = 0.0002822: For Case 1, total volume O 
= 0.04, with ply angle orientation (0°,,,,/ 90°/ 0°),,, (to scale) 


Of the two optimum solutions found for Case 1, the global optimum occurs 
when modes 1, 2, and 4 are simultaneous. This is not the expected result in which all 
four buckling modes occur simultaneously. The increase in size of the stiffener for the 
global optimum result allows the , (0°) layer to decrease down to the lower limit. This 
causes the stiffener to be relatively rig'd and the plate to be relatively weak in the loading 
direction. Therefore, for the global optimum case, a much greater portion of the load 
is carried by the stiffener. 
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Figure 6. Buckling Shapes for Case 1 with p = 0.0002822: For Case 1, total 


volume © = 0.04, with ply angle orientation of (0°,.,,,/ 90°/ 0°),,,. Re- 
presented above are the buckling shapes of the simultaneous overall 


plate buckling modes. 





ree 





2. Case 2: @ = 0.04 and Ply Orientation (0°,,,../ 0°/ 90°),. 
For this configuration, only one optimum solution is identified. It is listed in 
Table 2 below. The same results are presented graphically in Figures 7 and 8. This 
turns out to be the best possible configuration of those tested. Note that here all four 


modes buckle simultaneously. 


Table 2. RESULTS: @ = 0.04, PLY ANGLE (0°,¢,5/ 0°/ 90°) sy 


Buckling 
rs ee: 


0.00103 0.01385 0.06960 0.07350 | 0.0003949 








Figure 7. Geometry for Case 2 with p = 0.0003949: For Case 2, total volume © 
= 0.04, with ply angle orientation (0°,,,,,/ 0°/ 90°),,., (to scale) 
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Figure 8. Buckling Shapes for Case 2 with p = 0.0003949: For Case 2, total 
volume © = 0.04, with ply. angle orientation of (0°,,./0°/ 90°),,,. Re- 
presented above are the buckling shapes of the simultaneous overall 
plate buckling modes. 
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3. Case 3: © = 0.04 and Ply Orientation (90°, / 90°/ 0°),,. 
For this configuration, two different optimum solutions are identified. They are 


listed in Table 3 below. The same results are presented graphically in Figures 9-12. 


Table 3. RESULTS: @ = 0.04, PLY ANGLE (90° se44/ 90°/ O°) syy 


h Buckling 
Modes 
0.00100 0.01900 0.00100 0.00100 | 0.0000677 | 1 
0.00100 0.01538 0.02354 0.15323 | 0.0001804 









Figure 9. Geometry for Case 3 with p = 0.0000677: For Case 3, total volume © 
= 0.04, with ply angle orientation (90°,.,,,/ 90°/ 0°). 5» (to scale) 
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Buckling Shapes for Case 3 with p = 0.0000677: For Case 3, total 
volume © = 0.04, with py angle orientation of (90°,..../ 90°/ 0°),,.. 


Figure 10. 
Represented above are the buckling shapes of the simultaneous overall 


plate buckling modes. 














Figure 11. Geometry for Case 3 with p = 0.0001804: For Case 3, total volume 
© = 0.04, with ply angle onentation (90°,,,,,/ 90°/ 0°), (to scale) 


Of the two optima found, the global optimum occurs when three buckling 
modes occur simultaneously. No optimum was located that gave simultaneous buckling 
with all four modes. For this case, the stiffener fibers are in the 90° orientation so the 
stiffener adds much less strength to the plate than when the stiffener fibers are in the 0° 
orientation. Thus more of the total volume must be placed into the ¢, (0°) lamina. Be- 
cause of the volume of material in the /, lamina and the stiffener, there is an insufficient 
volume of material to allow the r, (90°, lamina to increase from the lower limit. Because 
of this, it is impossible to obtain simultaneous buckling with all four modes. 

For the local optima case where p = 0.0000677, the stiffened plate is trying to 
reduce itself down to a single lamina with 0° fiber orientation. From this configuration, 


any incremental increase in the stiffener size causes a decrease in buckling load. By ini- 
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tially increasing the stiffener size, the load share carried by the stiffener, p,, increases with 
the increase in the cross sectional area, s,. The initial increase in p, is greater than the 
increase in stiffness caused by the increase in s,, therefore, the buckling load initially 


decreases creating the local optimum. 
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Figure 12. Buckling Shapes for Case 3 with p = 0.0001804: For Case 3, total 
volume © = 0.04, with ply angle orientation of (90°,,,,,/ 90°/ 0°),,,. 
Represented above are the buckling shapes of the simultaneous overall 
plate buckling modes. 
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4. Case 4: © = 0.04 and Ply Orientation (90°... / 0°/ 90°), 
For this configuration, three different optimum solutions are identified. They 
are listed in Table 4 below. The same results are presented graphically in Figures 13-18. 


Table 4. RESULTS: @ = 0.04, PLY ANGLE (90° ,.,4/ 0°) 90°) syy 









Buckling 
Modes 


0.01900 | _0.00100__ | 0.00100_|_0.00100 Poeeor [1 
0.01636 0.00100 0.01845 0.14210 | 0.0001345 
0.00778 0.00816 0.02830 0.14310 | 0.0001702 








0.0000677: For Case 4, total volume 
(to scale) 


Figure 13. Geometry for Case 4 with p = 
© = 0.04, with ply angle orientation (90°,..,,/ 0°/ 90°),,. 
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Buckling Shapes for Case 4 with p = 0.0000677: For Case 4, total 
volume © = 0.04, with ply angle orientation of (90°,.,,,/ 0°/ 90°),,,: 
Represented above are the buckling shapes of the simultaneous overall 


Figure 14. 


plate buckling modes. 
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Figure 15. Geometry for Case 4 with p = 0.0001345: For Case 4, total volume 
© = 0.04, with ply angle orientation (90°,,,,,/ 0°/ 90°),,, (to scale) 


For this case, the global optimum occurs when all four buckling modes occur 
simultaneously. However, this global optimum is less than the global optimum of Case 
3. It appears as though the inner lamina (¢,) is the preferred location for the 0° fibers. 

The configuration which gives p = 0.0001345 corresponds closely to the Case 
3 global optimum. The only difference is that more 0° fiber orientation volume is re- 
quired when it is located at the outer layer (r,). The local optimum, p = 0.0000677, 
corresponds exactly to the Case 3 local optimum with the same load value. Here the 
same reasons apply for the existence of this local optimum. 


32 





Figure 16. 
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Buckling Shapes for Case 4 with p = 0.0001345: For Case 4, total 
volume © = 0.04, with ply angle orientation of (90°,./ 0°/ 90°) 
Represented above are the buckling shapes of the simultaneous overall 
plate buckling modes. 
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Figure 17. Geometry for Case 4 with p = 0.0001702: For Case 4, total volume 
© = 0.04, with ply angle orientation (90°,,,,,/0°/ 90°),,, (to scale) 








Figure 18. 





Buckling Shapes for Case 4 with p = 0.0001702: For Case 4, total 
volume © = 0.04, with ply angle orientation of (90°,,,,,/0°/ 90°),,.. 
Represented above are the buckling shapes of the simultaneous overall 
plate buckling modes. 
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5. Case 5: © = 0.015 and Ply Orientation (0°,.,,, / 90°/ 0°). 
For this configuration, two optimum solutions are identified. They are listed in 


Table 5 below. The same results are presented graphically in Figures 19-22. 


Table 5. RESULTS: 6 = 0.015, PLY ANGLE (0°,,,5,/ 90°/ 0°) yy 


Buckling 

sal 
0.00298 0.00411 0.01008 0.04029 
0.00576 0.00100 0.02720 0.02720 {| 0.0000124 camera 














Figure 19. Geometry for Case 5 with p = 0. 0000156: For Case 5, total volume 


© = 0.015, with ply angle orientation (0°,,,,,/ 90°/ 0°)... (to scale) 
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Figure 20. 





Buckling Shapes for Case 5 with p = 0.0000156: For Case 5, total 

volume © = 0.015, with-ply angle onentation of (0°,,,,./ 90°/ 0°),,,. 
Represented above are the buckling shapes of the simultaneous overall 
plate buckling modes. 
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Figure 21. Geometry for Case 5 with p = 0.0000124: For Case 5, total volume 
© = 0.015, with ply angle orientation (0°,.20/ 90°/ 0°). (to scale) 


Here the global optimum is found when the four buckling modes occur sircul- 
taneously. This global optimum corresponds to the global optimum found in Case 1. 
This follows since Case 5 is the same as Case 1, except for the smaller volume of material 
available. The local optimum found with p = 0.0000124 is not a uniquely determined 
solution. The total stiffener voiume is uniquely defined but not the individual dimen- 
sions of the stiffener height and width. This occurs because the first buckling mode (the 
only active mode in this case) does not cause any flexure of the stiffener. The active 
buckling shape only twists the stiffener so that only its volume and not its cross sectional 
shape is important to the overall buckling stiffness. 
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Figure 22. Buckling Shapes for Case § with p = 0.0000124: For Case 5, total 
volume © = 0.015, with ply angle orientation of (0°,../ 90°/ 0°),.»- 


Represented above are the buckling shapes of the simultaneous overall 


plate buckling modes. 
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6. Case 6: @ = 0.015 and rly Orientation (0°... / 0°/ 90°), , 
For this configuration, only one optimum solution is identified. It is listed ir 
Table 6 below. The same results are presented graphically in Figures 23-24. This case 
is simular to Case 2 except here there is insufficient material available to produce simul- 


taneous buckling of all four modes. 


Table 6. RESULTS: 6 = 0.015, PLY ANGLE (0° 5.44, / 0°/ 90°) yy 


Buckling 
Modes 


0.00100 0.00577 0.01814 0.03986 1 











Figure 23. Geometry for Case 6 with p = 0.0000179: For Case 6, total volume 


© = 0.015, with ply angle orientation (0°,,,,./ 0°/ 90°), (to scale) 
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Figure 24. Buckling Shapes for Case 6 with p = 0.0000179: For Case 6, total 
volume © = 0.015, with ply angle orientation of (0°,.,,/ 0°/ 90”) ym’ 
Represented above are the buckling shapes of the simultaneous overall 
plate buckling modes. 
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7. Case 7: @ = 0.015 and Ply Orientation (90°... / 90°/ 0°),,. 
For this configuration, two different optimum solutions are identified. They are 
listed in Table 7 below. The same results are presented graphically in Figures 25-28. 





Table 7. RESULTS: @ = 0.015, PLY ANGLE (90° c.4/ 90°/ 0°) say 
Buckling 


Sa ies |e 


0.00100 0.60650 0.00100 0.00100 | 00000036 | 1 | 
0.00100 0.00584 0.00669 0.09761 | 0.0000130 










Figure 25. Geometry for Case 7 with p = 0.0000036: For Case 7, total volume 
© = 0.015, with ply angle oricntation (90°,,,,,/ 90°/ 0°), 5» (to scale) 
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Buckling Shapes for Case 7 with p = 0.0000036: For Case 7, total 
volume © = 0.015, with ply angle orientation of (90°,,,,,/ 90°/ 0°),,.- 


Figure 26. 
Represented above are the buckling shapes of the simultaneous overall 


plate buckling modes. 
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Figure 27. Geometry for Case 7 with p = 0.0000130: For Case 7, total volume 
© = 0.015, with ply angle orientation (90°,,,,,/ 90°/ 0°),,,, (to scale) 


The optima of Case 7 are closely related to the optima of Case 3. The only 
difference in the two cases is the total volume of material available. Here the same rel- 
ative magnitude between the two different optima are present with the global optima 
again being the configuration with three buckling modes occurring simultaneously. As 
in Case 3, here too, there is insufficient material available to obtain the case where all 
four buckling modes occur simultaneously. 
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Figure 28. Buckling Shapes for Case 7 with p = 0.0000130: For Case 7, total 
volume © = 0.015, with ply angle orientation of (90°,.,../ 90°/ 0°),,,.. 
Represented above are the buckling shapes of the simultaneous overall 


plate buckling modes. 
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8. Case 8: © = 0.015 and Ply Orientation (90°... / 0°/ 90°), 
For this configuration, two different optimum solutions are identified. They are 
listed in Table 8 below. The same results are presented graphically in Figures 29-32. 


Table 8. RESULTS: 6 = 0.015, PLY ANGLE (90°,,,5,/ 0°/ 90°), 


Buckling 

Modes 
0.00650 0.00100 0.00100 0.00100 | 0.0000036 | 1 | 
0.00329 0.00345 0.00802 0.09386 NCTE 









Figure 29. Geometry for Case 8 with p = 0.0000036: For Case 8, total volume 
© = 0.015, with ply angle orientation (90°,,,,,/ 0°/ 90°),,., (to scale) 
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Figure 30. Buckling Shapes for Case 8 with p = 0.0000036: For Case 8, total 
volume © = 0.015, with ply angle orientation of (90°,..,/ 0°/ 90") .,5: 
Represented above are the buckling shapes of the simultaneous overall 
plate buckling modes. 


47 





Pte 








Figure 31. Geometry for Case 8 with p = 0.0000125: For Case 8, total volume 
© = 0.015, with ply angle orientation (90°,,,,,/ 0°/ 90°), (to scale) 


As in Case 4, the global optimum solution for Case & occurs when all four 
buckling modes occur simultaneously. Also, there is a similar |. .al optimum to that of 
Case 4, which has only one active buckling mode. This case is present for the same 
reasons as given in Case 3. One difference between Case 8 and Case 4 is that no local 
optimum in which three buckling modes occur simultaneously could be found for Case 
8 as was found in Case 4. 





Figure 32. 
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Buckling Shapes for Case 8 with p = 0.0000125: For Case 8, total 
volume © = 0.015, with ply angle orientation of (90°, / 0°/ 90°),,.. 


Represented above are the buckling shapes of the simultaneous overall 


plate buckling modes. 
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B. OPTIMIZATION EFFICIENCY 


To compare the efficiencies of the three optimization routines, two separate criteria 


are used. These criteria are: the total execution time required to converge to an opti- 
mum solution; and the ability of an optimizer to converge to a solution when the initial 
guess is far from the solution (i.e. globally convergent). In order to compare the relative 
efficiencies of the three optimization methods, the exact same problem is solved with 
each of the optimization schemes. For each run, the same initial guess is used and the 
actual execution time is measured. These results are presented in Table 9 below. 


Table 9. RUN TIMES: © = 0.04, PLY ANGLE (0° ,45,/ 90°/ 0°) yy 


Initial Guess 
pif | 6 | 4 | p _| ae 
cution 




























Opti- 

mization 7.00100 | 0.10000 [0.10060 | 0:6001000 | Time in 

oe Te oe 
seconds 


ip aS eee td 
SS ET TT TE 






From these results, 1t can be seen that for this specific design problem, the IMSL 
routine is clearly the most efficient in terms of the execution time. The relative speed 
of convergence for the IMSL routine is due to the high degree of nonlinearity in the 
design problem. Because the IMSL routine uses sequential quadratic programming as 
its search algorithm, which is based on a second order derivative, it allows each iteration 
in the optimum search to be more profitable than those of the DOT algonthms. Both 
the DOT methods (MFD and SLP) have search directions that are based on first order 
derivatives. Although the search direction calculation for both DOT methods is simpler, 
it is not as profitable in each iteration in the optimum search as that of the IMSL 
method. Thus, more iterations are required for optimum solution convergence with the 
DOT methods. The savings made in the calculation time for each search direction in the 
DOT code is not nearly enough to make up for the additional calculation time rcuuired 
for the extra iterations in the optimization Therefore, due to the highly nonlin. na- 
ture of this problem, the IMSL routine is 2 efficient in terms of execution tin 


50 





lene 














The second measure of the optimizer’s efficiency is its ability to be globally conver- 
gent. That is, to converge to some solution from almost any starting point [Ref. 32: p. 
5]. Since the location and number of optima are unknown at the start of the optimiza- 
tion process, the initial guesses are scattered throughout the design space in an attempt 
to identify all possible optima. Because the IMSL routine is based on a second order 
search, it is better suited to converge when the initial guess is relatively far away from 
an optuma. Both the DOT methods were unable to converge to an optima on many runs 
when the initial guess was relatively far from any optima. Qualitatively, the IMSL rou- 
tine produced actual optimum convergence ten fold more often than either of the DOT 
methods. Therefore, from the standpoint of globally convergent efficiency, the IMSL 
routine is clearly superior for this specific nonlinear optimization problem. 
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Vv. CONCLUSIONS 


Many times it is difficult to draw general conclusions for design problems. This case 
is no exception, however, there are three main areas where specific conclusions for this 
particular design problem can be drawn. These areas are: the general character of the 
design problem, the best configuration of the investigated stacking sequences; and the 
relative efficiencies of the three optimization routines used. 

The first general conclusion drawn is the character of the design problem for the 
maximum buckling load of the stiffened composite plate. The design optimization 
problem turns out to be a highly nonlinear problem. For almost all cases investigated, 
there were multiple optimum solutions indicating a nonlinear solution to the design op- 
timization problem. Some cases had as many as three different optima. Even the case 
where only one optimum solution is identified cannot be assumed to be an indication 
of a linear problem. This must be viewed as the only solution identified and not as the 
only solution that exists. Additionally, these solutions were all located in a limited de- 
sign space. One would expect even more optima to be found by increasing the size of 
the design space. Finally, the order of the governing partial differential equation would 
lead one to assume a highly nonlinear problem without examining the results. The re- 
sults merely confirm this assumption. 

The second area of general conclusions is that of the best design configuration. For 
both total volume constraints (QO = 0.04 and 0.015) examined, the best stacking config- 
uration was the (0°,.2,/0°/ 90°),,, used in cases 2 and 6. The use of the 0°,,,,, stiffener 
fiber orientation was always superior to the 90°,,,, stiffener fiber orientation regardless 
of the laminae ply orientation. This led to the conclusion that the 0°,,,, stiffener fiber 
orientation is the most important factor in the design optimization. The laminae fiber 
orientation of (0°/90°),,, was best when used with the 90°,,,, stiffener fiber orientation. 
However, when used with the 0°,,,, stiffener fiber orientation, the (90°/0°),,, laminae fi- 
ber orientation was best. This indicates that for best results, the inner lamina (¢, layer) 
should have its fiber orientation perpendicular to that of the stiffener fiber orientation. 

The third area of general conclusions is the relative efficiencies of the three opti- 
mization routines used. Because of the highly nonlinear nature of this specific design 
optimization problem, the IMSL subroutine DNCONG proved to be the most efficient, 


both in terms of convergence speed and in its ability to be globally convergent (i.e. 


52 


rlave 


converge to an optimum from almost any starting point). Because the IMSL routine is 
based on a second order search direction, it is better suited to optimize the highly non- 
linear design problem. In contrast to this, both of the DOT methods have search di- 
rections that are based on first order search directions that make them ill-suited for such 
a highly nonlinear optimization problem. 

It is important to note that the conclusions of this research are specific to the spe- 
cific design optimization problem investigated. Generalizations of these results to other 
problems must be examined closely due to the nonlinear nature of the problem. Any 
design optimization problem of a nonlinear system must be closely scrutinized by the 
designer to ensure that all important aspects of the design are considered. Computer- 
aided design optimization can be a very powerful tool for a designer, however, caution 
must be exercised to ensure that a local optimum design is not selected as the best pos- 


sible design without investigating the entire design space for other optima. 
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APPENDIX A. DESIGN ANALYSIS SUBROUTINES 


SUBROUTINE EIGANL4 


ae ce a ee me ee a cs ee a ee a a a nw ee ee ee oe a ee ee 
SSS S255 5552552555 5525552555555 5552555552552 5552555225552 2>>=S=>== 


LINEAR FINITE ELEMENT ANALYSIS OF LAMINATED COMPOSITE PLATE 
(16 DOF RECTANGULAR PLATE ELEMENT) 
3 LOWEST EIGENVLAUES AND THEIR DERIVATIVES WRT DESIGN VARIABLES 


re me ee ee ee a ee a ee oe ee we ew a a we ee a a a ee 
SSS SSS S522 5525555255205 5552525555255 >5ES 52S =5=S5=SS5S—> 


IMPLICIT REAL*8(A-H,0~Z) 

COMMON/ELEMNT/ EK(16,16,6), EKT( 16,16), EG(16,16) 

COMMON/ELEMNB/ ARSTF(4,4), ARGEM(4,4) 

COMMON/PROPT1/ E£1(15), E2(15), G12(15), PO12(15), ANGL(15), 2(16) 
» NOD(150,4), NX, NY, NSDF, LSDF(150) 

COMMON/PROPT2/ DM(3,3) 

COMMON/PROPT3/ EBEAM, ABEAM, BMINER, EPLT, APLT, THICK, NODB(18,2) 

COMMON/PROPT4/ HT(15), AAA, BBB, NLYR, NLYR1, NLYR2, NLYR3 

COMMON/TOTAL / GK(150,150) ,GG( 150,150) : 

COMMON/BAND/ BANDK(45,150), BANDKR(45,150), BANDG(45,150) 

COMMON/EIG/NBAND, NREQ, NEQ 

COMMON/EIG2/ EGVC(150,5), EGVCO(150,5), EGVL(5,4), EGVLO(5,4), 
DTHK(50,5), DGD(5), EGVCB(150,5), IDMOD 

DIMENSION DBANDK(45,150), WK(150), WK1( 150), WK2(150) 

COMMON/BEAM/ E1B, E2B, G12B, PO12B, BBM, HBM 

COMMON/BEAM2/ DBTHK(50), PBEAM 

COMMON/NORMAL/ VNORMAL 


FREE 


CALL PROPTY 
CALL ELESTF 
CALL ASSEMB 
CALL EIGEN 
IDMOD = 1 


DO 11 I=1,3 
EGVL(I,1) = EGVL(1I,1)/VNORMAL 
CONTINUE 


SAVE K, G 


NBAND1 = NBAND + 1 
DO 150 I=1, NREQ 
DO 140 J=1, NBAND1 
DBANDK(J,1) = BANDKR(J,I) - EGVL(1,1)*BANDG(J,I) 
CONTINUE 
CONTINUE 




















C === OBTAIN DERIVATIVES OF P WRT T 
C 


DO 1305 J=1, NREQ 
WK(J) = EGVC(J,1) 
1305 CONTINUE 
Lc = 45 
CALL BNDMUL(BANDG ,WK ,NBAND1,NREQ,WK2,LC) 
DGD(1) = 0.Dd0 
DO 1320 J=1, NREQ 
DGD(1) = DGD(1) + WK(J)*WK2(J) 
1320 CONTINUE 
IFC IDMOD.GE.1) THEN 
DO 1307 J=1, NREQ 
WK1(J) = EGVC(J,2) 
1307 CONTINUE 
LC = 45 
CALL BNDMUL( BANDG ,WK1,NBAND1,NREQ ,WK2,LC) 
DGD(2) = 0.D0 
DO 2320 J=1, NREQ 
DGD(2) = DGD(2) + WK1(J)*WK2(J) 


2320 CONTINUE 

C 

C xk FOR 3RD EIGENVECTOR 
Cc 


CALL BNDMUL(BANDG ,EGVC(1,3) ,NBAND1,NREQ,WK2,LC) 
DGD(3) = 0.D0 
DO 4320 J=1, NREQ 

DGD(3) = DGD(3) + EGVC(J,3)*WK2(J) 


4320 CONTINUE 
END IF 

Cc 

Cc 

C 


DEL = 1.0D-05 
DO 380 IV=1, NLYR2 
IF(IV.EQ.NLYR1) THEN 
BBM = BBM + DEL 
ELSE IF(IV.EQ.NLYR2) THEN 
HBM = HBM + DEL 
ELSE 
HT(IV) = HT(IV) + DEL 
END IF 
CALL PROPTY 
CALL ELESTF 
CALL ASSEMB 


C 
Cres GET DK 
Cc 
DO 300 J=1, NREQ 
DO 280 K=1, NBAND1 ; 
BANDKR(K,J) = (BANDKR(K,J) - EGVL(1,1)*BANDG(K,J) 

i - DBANDK(K,J))/DEL 
280 CONTINUE 
300 CONTINUE 
C 
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C --- 


3320 


aaAanNnaANa 


GET 


U*DK*U 


LC = 45 

CALL BNDMUL(BANDKR ,WK ,NBAND1,NREQ,WK2,LC) 
DTHK(IV,1) = 0.DO 

DO 320 J=1, NREQ 


DTHK(IV,1) = DTHK(IV,1) + WK(J)*WK2(J) 


CONTINUE 
DTHK(IV,1) = DTHK(IV,1)/DGD(1) 


GET 


U*DK*U 


IF( IDMOD.GE.1) THEN 


END 


LC = 45 
CALL BNDMUL( BANDKR ,WK1,NBAND1,NREQ ,WK2,LC) 
DTHK(IV,2) = 0.D0 
DO 3320 J=1, NREQ 

DTHK(IV,2) = DTHK(IV,2) + WK1(J)*WK2(J) 
CONTINUE 
DTHK(IV,2) = DTHK(IV,2)/DGD(2) 


CALL BNDMUL(BANDKR ,EGVC(1,3) ,NBAND1 ,NREQ ,WK2,LC) 
DTHK(IV,3) = 0.D0 
DO 5320 J=1, NREQ 

DTHK(IV,3) = DTHK(IV,3) + EGVC(J,3)*WK2(J) 
CONTINUE 
DTHK(IV,3) = DTHK(IV,3)/DGD(3) 
IF 


ripe 


ORIGINAL HT 


IF(IV.EQ.NLYR1) THEN 


BBM = BBM - DEL 


ELSE IF(IV.EQ.NLYR2) THEN 


HBM = HBM - DEL 


ELSE 


END 


HT(1IV) = HT(IV) - DEL 


IF 


CONTINUE 


DO 4694 I=1,NLYR2 
DO 4694 J=1,3 


DTHK(1,J) = VNORMAL*DTHK‘I,J) 
CONTINUE 


DO 1111 I=1,3 
EGVL(I,1) = EGVL(1I,1)*VNORMAL 


RETURN 
END 


CONTINUE 
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SUBROUTINE INPUTB 


C 
Cc 
IMPLICIT REAL*8(A-H,0-Z) 

COMMON/PROPTO/ THETA, ALEN, BLEN 

COMMON/PROPT1/ E£1(15), E2(15), G12(15), PO12(15), ANGL(15), 2(16) 

2 »NOD(150,4), NX, NY, NSDF, LSDF(150) 

COMMON/PROPT4/ HT(15), AAA, BBB, NLYR, NLYR1, NLYR2, NLYR3 

COMMON/BEAM/ E1B, E2B, G12B, PO12B, BBM, HBM 
c 
C READ IN GEOMETRIC PARAMETERS AND PHYSICAL PROPERTIES 
Cc 


READ(5,*) ALEN, BLEN 

READ(5,*) NLYR,THETA 

NLYR1 = NLYR + 1 

NLYR2 = NLYR + 2 

NLYR3 = NLYR + 3 

DO 150 I=1, NLYR 

READ(S,*) E1l(1I), E2(1I), G12(1), PO12(I), ANGL(I) 

150 CONTINUE 

READ(5,*) NX, NY 

READ(5,*) E1B, E2B, G12B, PO12B 

RETURN 


rare 





SUBROUTINE INPUTC 





qa aang 


IMPLICIT REAL*8(A-H,0-Z) 

EXTERNAL GRMESH 

COMMON/PROPTO/ THETA, ALEN, BLEN 

COMMON/PROPT1/ E1(15), E2(15), G12(15), PO12(15), ANGL(15), Z(16) 
* » NOD(150,4), NX, NY, NSDF, LSDF(150) 

COMMON/PROPT3/ EBEAM, ABEAM, BMINER, EPLT, APLT, THICK, NODB(18,2) 
COMMON/PROPT4/ HT(15), AAA, BBB, NLYR, NLYR1, NLYR2, NLYR3 
COMMON/EIG/NBAND, NREQ, NEQ 

COMMON/EIG2/ EGVC(150,5), EGVCO(150,5), EGVL(5,4), EGVLO(5,4), 
* DTHK(50,5), DGD(5), EGVCB(150,5), IDMOD 
COMMON/BEAM/ E1B, E2B, G12B, PO12B, BBM, HBM 


Z(1) = 0.D0 
DO 150 I=1, NLYR 
ZCI+1) = ZCI) + HTCI) 
150 CONTINUE 
THICK = Z(NLYR1)*2.0D0 
APLT = BLEN*THICK 


NEL = NX*NY 

CALL GRMESH(NNM ,NEM) 
C **% BEAM CONNECTIVITY MATRIX *** 
Cc 

NX1 = NX + 1 

NY1 = NY + 1 
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NODO = NY/2*NX1 + 1 
DO 124 I=1, NX 
NODN = NODO + 1 
NODB(I,1) = NODO 
NODB(I,2) = NODN 
NODO = NODN 
124 CONTINUE 


C 
C *** BOUNDARY CONDITIONS OF SIMPLE SUPPORT *** 
c 

NSDF = 12 + 4*(NX-1) + 4*(NY-1) 

N = 0 

DO 400 K=1, NY1 

DO 300 J=1, NX1 
IF(K.NE.1.AND.K.NE.NY1) THEN 
IF(J.NE.1.AND.J.NE.NX1) THEN 
GO TO 300 


4*((K-1)*NX1 + J - 1) + 1 
200 CONTINUE 
END IF 
ELSE 
IF(J.NE.1.AND.J.NE.NX1) THEN 
DO 250 I=1, 3, 2 


N=N+1 
LSDF(N) = 4*((K-1)*NX1 + J - 1) +I 
250 CONTINUE 
ELSE 
DO 270 I=1, 3 
N=Nt1 
LSDF(N) = 4*((K-1)*NX1 + J - 1) +1 
270 CONTINUE 
END IF 
END IF 
300 CONTINUE 
400 CONT .NUE 
C 
AAA = ALEN/NX 
BBB = BLEN/NY 
NEQ = 4*(NX+1)*(NY41) 
NREQ = NEQ - NSDF 
C 
C --- CLEAR EGVCB 
C 
DO 722 I=1, 3 
DO 721 J=1, NREQ 
EGVCB(J,1) = 0.D0 
721 CONTINUE 
722 CONTINUE 
RETURN 
Cc 
C 
Cc 
C 
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SUBROUTINE PROPTY 


IMPLICIT REAL*8(A-H,0-Z) 
COMMON/PROPT4/ HT(15), AAA, BBB, NLYR, NLYR1, NLYR2, NLYR3 
COMMON/PROPT1/ E1(15), E2(15), Gi2(15), PO12(15), ANGL(15), Z(16) 


* » NOD(150,4), NX, NY, NSDF, LSDF(150) 


COMMON/PROPT2/ DM(3,3) 

COMMON/PROPT3/ EBEAM, ABEAM, BMINER, EPLT, APLT, THICK, NODB(18,2) 
COMMON/BEAM/ E1B, E2B, G12B, PO12B, BBM, HBM 

DIMENSION DANGL(15) 


Z(1) = 0.D0 
DO 5 I=1, NLYR 

Z2(I+1) = Z2(1I) + HT(I) 
CONTINUE 


DO 10 I=1, 3 
DO 10 J=1, 3 

DM(I,J)=0.0 
CONTINUE 
AM11 
AM12 
AM13 
AM22 
AM23 
AM33 


.DO 
.DO 
.DO 
.DO 
.DO 
-DO 


ecrte 


ooooo0o°o 


DO-LOOP FOR NUMBER OF LAYERS 
DO 30 I=1,NLYR 
CALCULATION OF THE REDUCED STIFFNESSES 


PO21=E2(1)*P012(1)/E1(1) 
BPO=1.0-P012(1)*P021 
Q11=Ei(1I)/BPO 
Q12=P012(1)*E2(1)/BPO 
Q22=E2(1)/BPO 
Q66=G12(T) 


CALCULATION OF THE TRANSFORMED REDUCED STIFFNESSES 


DANGL(I)=3.141592D0*ANGL(I)/180.D0 
SC=DSIN( DANGL( I) )*DCOS(DANGL(I)) 
S2=(DSIN(DANGL( TI) ) )**2 
C2=(DCOS(DANGL( 1) ) )**2 


QB11=Q11*C2*C2+2 .0*(Q12+2 . 0*Q66 ) *SC*SC+Q22*S2*S2 
QB22=Q11*S2*S2+2 .0*(Q12+2.0*Q66 )*SC*SC+Q22*C2*C2 
QB12=(Q11+Q22-4. 0*Q66)*SC*SC+Q12*(S2*S2+C2*C2) 
QB16=(Q11-Q12-2 .0*Q66)*SC*C2+(Q12-Q22+2 . 0*Q66)*SC*S2 
QB26=(Q11-Q12-2.. 0*Q66 )*SC*S2+(Q12-Q22+2 . O*Q66) *SC*C2 
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30 


QB66=(Q114+Q22-2 .0*Q12-2.. 0*Q66) *SC*SC+Q66*(S2*S2+C2*C2) 
CALCULATION OF THE EXTENSIONAL STIFFNESSES 


AM11=QB11*((Z(1I+1))-(Z(1)))}*2.0D0 + AM11 
AM12=QB12*((Z(I+1))-(Z2(1)))*2.0D0 + AM12 
AM22=QB22*((Z(I+1))-(Z(1)))*2.0D0 + AM22 
AM13=QB16*((Z(I+1))-(2Z(1)))*2.0D0 + AM13 
AM33=QB66*((Z(I+1))-(Z(1)))*2.0D0 + AM33 
AM23=QB26*( (ZC 1+1))-(Z(1)))*2.0D0 + AM23 


CALCULATION OF THE BENDING STIFFNESSES 


DM(1,1)=QB11*((Z( I+1) )**3-(Z(1) )**3)/1.5D0 
DM(1,2)=QB12*( (ZC 1+1) )**3-(Z(1) )**3)/1.5D0 
DM(2,2)=QB22*( (ZC I+1) )**3-(Z2(1) )**3)/1.5D0 
DM(1,3)=QB16*((Z(I+1) )**3-(Z(1))**3)/1.5D0 
DM(3,3)=QB66*((Z(I+1) )**3-(Z(1I) )**3)/1.5D0 
DM(2,3)=QB26*( (ZC I+1) )**3-(Z(1) )**3)/1.5D0 

CONTINUE 

AM21 = AM12 

AM31 = AM13 

AM32 = AM23 


DM(1,1) 
DM(1,2) 
DM(2,2) 
DM(1,3) 
DM(3,3) 
DM(2,3) 


tee eet 


FIND E OF THE PLATE *** 


wegen 


BLEN = BBB*NY 
THICK = Z(NLYR1)*2.0D0 
APLT = BLEN*THICK 


DENO = AM11*AM22*AM33 + AM12*AM23*AM31*2.D0 - AM13**2*AM22 - 
AM12**2*AM33 - AM11*AM23**2 

EPLT = DENO/(AM22*AM33 - AM23**2)/THICK 

DM(2,1) = DM(1,2) 


DM(3,1) = DM(1,3) 

DM(3,2) = DM(2,3) 

ABEAM = 2.D0*BBM*HBM 

EBEAM = E1B 

THIK2 = THICK/2.D0 

BMINER = 2.D0*BBM/3.D0*( (THIK2+HBM)**3-THIK2**3) 
RETURN 

END 


SUBROUTINE ELESTF 


IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION GAUSS4(4), GAUSS3(3), WT4(4), WT3(3), BG(16) 
COMMON/ELEMNT/ EK(16,16,6), EKT(16,16), EG(16,16) 

COMMON/ELEMNB/ ARSTF(4,4), ARGEM(4,4) 

COMMON/PROPT2/ DM(3,3) 

COMMON/PROPT3/ EBEAM, ABEAM, BMINER, EPLT, APLT, T’ CK, NODB(18,2) 
COMMON/PROPT4/ HT(15), AAA, BBB, NLYR, NLYR1, NLYF NLYR3 
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DIMENSION BK(16,3), WK1616(16,16), WK163(16,3), WK316(3,16), 
* WK33(3,3), DV(6) 


DATA GAUSS4/ .0694318D0 ,0. 3300095D0 ,0.6699905D0 ,0.9305681D0/ 
DATA GAUSS3/.1127017D0,0.5D0,0.88729833D0/ 

DATA WT4/0.1739274D0 ,0.3260725D0 ,0.3260725D0 ,0.1739274D0/ 
DATA WT3/0.2777777D0,0 .4444444D0 ,0.2777777D0/ 


INTERPOLATION FUNCTIONS AND THEIR DERIVATIVES *** 


FNI(P) = 1.D0 - 3.DO0*P*P + 2.D0*P*P*P 
FN2(P) = (1.D0 - 2.D0*P + P*P)*P 
FN3(P) = (3.D0 - 2.D0*P)*P*P 
FN4(P) = P*P*(1.D0 - P) 

DFNI(P) = P*(-6.D0 + 6.D0*P) 
DFN2(P) = 1.D0 - 4.D0*P + 3.D0*P*P 
DFN3(P) = 6.DO0*P*(1.D0 - P) 
DFN4(P) = P*(2.D0 - 3.DO*P) 
DDFN1(P) = -6.D0 + 12.D0*P 
DDFN2(P) = -4.D0 + 6.DO*P 

DDFN3(P) = 6.D0 - 12.DO0*P 

DDFN4(P) = 2.D0 - 6.DO*P 


COMPONENT OF D MATRIX 


eegit 


DV(1) = DM(1,1) 
DV(2) = DM(1,2) 
DV(3) = DM(1,3) 
DV(4) = DM(2,2) 
DV(5) = DM(2,3) 
DV(6) = DM(3,3) 


INITIALIZE THE ELEMENT MATRICES AND VECTORS 


DO 20 J=1, 16 
DO 20 I=1, 16 
EG(I,J) = 0.D0 
EXT(I,J) = 0.D0 
CONTINUE 


DO 30 K=1, 6 

DO 30 J=1, 16 
I=1, 16 
I 


»J,K) = 0.D0 


DO 30 
EK( 
CONTINUE 


NUMERICAL INTEGRATION BY GAUSS QUADRATURE 
DO 500 NI=1,3 
XI=GAUSS3(NI) 


DXF1 = DFNi(XI) 
DXF2 = DFN2(XI) 


61 








250 


DXF3 = DFN3(XI) 
DXF4 = DFN4(XI) 
DO 500 NJ=1,4 
ETA=GAUSS4(NJ) 
EF1 = FNI(ETA) 
FN2(ETA) 
FN3(ETA) 
FN4(ETA) 


is] 
i> | 
w 
wn nu 


DXF 1*EF1/AAA 

DXF 1*EF2*BBB/AAA 
-DXF2*EF1 
DXF2*EF2*BBB 


ow 
Q 
py 
to 
— 
nouw t 


DXF3*EF1/AAA 
DXF3*EF2*BBB/AAA 
DXF4*EF1 
-DXF4*EF2*BBB 


w 
Q 
~ 
n 
Ne 
uuu tb 


BG(9) = DXF3*EF3/AAA 
BG(10) = -DXF3*EF4*BBB/AAA 
BG(11) = DXF4*EF3 

BG(12) = DXF4*EF4*BBB 


BG(13) 
BG( 14) 
BG(15) 
BG( 16) 


DXF 1*EF3/AAA 
-DXF1*EF4*BBB/AAA 
-DXF2*EF3 
-DXF2*EF4*BBB 


DO 250 J=1, 16 
DO 250 I[=1, 16 
EG(I,J) = EG(I,J) + WT3(NI)*WT4(NJ)*BG(I)*BG(J) 
CONTINUE 


500 CONTINUE 


Cc 
C *** MULTIPLY THE LOAD FACTOR *** 


C 


555 


FACTOR = EPLT*THICK/(EPLT*APLT + EBEAM*ABEAM)*AAA*BBB 
DO 555 J=1, 16 
DO 555 I=1, 16 
EG(I,J) = EG(I,J)*FACTOR 
CONTINUE 


DO 600 NI=1, 4 
XI=GAUSS4(NI) 
XF1 = FN1(XI) 
XF2 = FN2(XI) 
XF3 = FN3(XI) 
XF4 = FN4(XI) 
DXF1 = DFN1(XI) 
DXF2 = DFN2(XI) 
DXF3 = DFN3(XI) 
DXF4 = DFN4(XI) 
DDXF1 = DDFN1(XI) 


huun 


ne 
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reg't 


DDXF2 
DDXF3 
DDXF4 


DDFN2(XI) 
DDFN3(XI) 
DDFN4(XI) 


DO 600 NJ=1, 4 
ETA=GAUSS4(NJ) 


EF1 
EF2 
EF3 
EF4 
DEF 1 
DEF2 
DEF3 
DEF4 
DDEF1 
' DDEF2 
DDEF3 
DDEF4 


Te | 


BK(1,1) 
BK(2,1) 
BK(3,1) 
BK(4,1) 


BK(5S,1) 
BK(7,1) 
BK(8,1) 


BK(9,1) 

BK(10,1) 
BK(11,1) 
BK(12,1) 


BK(13,1) 
BK(14,1) 
BK(15,1) 
BK( 16,1) 


BK( 1,2) 
BK( 2,2) 
BK( 3,2) 
BK( 4,2) 


BK( 5,2) 
BK( 6,2) 
" BK( 7,2) 
BK( 8,2) 


BK( 9,2) 
BK(10,2) 
BK(11,2) 
BK( 12,2) 


BK( 13,2) 
BK(14,2) 
BK(15,2) 


FN1(ETA) 
FN2(ETA) 
FN3(ETA) 
FN4(ETA) 
DFN1(ETA) 
DFN2(ETA) 
DFN3(ETA) 
DFN4(ETA) 


DDFN1(ETA) 
DDFN2(ETA) 
DDFN3(ETA) 
DDFN4(ETA) 


uuu wt 


DDXF 1*EF1/AAA/AAA 
DDXF 1*EF2*BBB/AAA/ AAA 
-DDXF2*EF1/AAA 

DDXF 2*EF2*BBB/AAA 


DDXF3*EF1/AAA/AAA 
DDXF3*EF2*BBB/AAA/AAA 
DDXF4*EF1/AAA 
~DDXF4*EF2*BBB/AAA 


DDXF3*EF3/AAA/AAA 
-DDXF3*EF4*BBB/AAA/AAA 
DDXF4*EF3/AAA 
DDXF4*EF4*BBB/AAA 


DDXF1*EF3/AAA/AAA 
-DDXF1*EF4*BBB/AAA/AAA 
-DDXF2*EF3/AAA 
-DDXF2*EF4*BBB/AAA 


XF 1*DDEF1/BBB/BBB 
XF1*DDEF2/BBB 
-XF2*DDEF 1*AAA/BBB/BBB 
XF2*DDEF2*AAA/ BBB 


XF3*DDEF 1/BBB/BBB 
XF3*DDEF2/BBB 
XF4*DDEF1*AAA/BBB/BBB 
-~XF4*DDEF2*AAA/ BBB 


XF3*DDEF3/3BB/BBB 
~XF3*DDEF4/ BBB 
XF4*DDEF3*AAA/BBB/BBB 
XF4*DDEF4*AAA/BBB 


XF1*DDEF3/BBB/BBB 


-XF1*DDEF4/BBB 
-XF2*DDEF3*AAA/BBB/BBB 
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BK(16,2) = -XF2*DDEF4*AAA/BBB 


G 
BK( 1,3) = 2.DO*DXF1*DEF1/AAA/BBB 
BK( 2,3) = 2.DO*DXF1*DEF2/AAA 
BK( 3,3) = -2.D0*DXF2*DEF1/BBB 
BK( 4,3) = 2.DO*DXF2*DEF2 
c 
BK( 5,3) = 2.DO*DXF3*DEF1/AAA/BBB 
BK( 6,3) = 2.D0*DXF3*DEF2/AAA 
BK( 7,3) = 2.DO*DXF4*DEF1/BBB 
BK( 8,3) = -2.DO*DXF4*DEF2 
c 
BK( 9,3) = 2.D0*DXF3*DEF3/AAA/BBB 
BK(10,3) = -2.DO*DXF3*DEF4/AAA 
BK(11,3) = 2.DO*DXF4*DEF3/BBB 
BK(12,3) = 2.DO*DXF4*DEF4 
c 
BK(13,3) = 2.D0*DXF1*DEF3/AAA/BBB 
BK(14,3) = -2.DO*DXF1*DEF4/AAA : 
BK(15,3) = -2.DO*DXF2*DEF3/BBB 
BK(16,3) = -2.D0*DXF2*DEF4 
c 
C wie GET TRANSPOSE OF BK 
c 
DO 570 I=1, 3 
DO 570 J=1, 16 : 
WK316(I,J) = BK(J,1) : 
570 CONTINUE 
c 


DO 580 KK=1, 6 

DO 560 JL=1, 3 

DO 560 IL=1, 3 
WK33(IL,JL) = 0.D0 
560 CONTINUE 

IF(KK.EQ.1) THEN 
WK33(1,1) = 1.D0 

ELSE IF(KK.EQ.2) THEN 


WK33(1,2) = 1.D0 
WK33(2,1) = 1.D0 

ELSE IF(KK.EQ.3) THEN 
WK33(1,3) = 1.D0 
WK33(3,1) = 1.D0 

ELSE IF(KK.EQ.4) THEN : 
WK33(2,2) = 1.D0 

ELSE IF(KK.EQ.5) THEN 
WK33(2,3) = 1.D0 : 
WK33(3,2) = 1.D0 

ELSE IF(KK.EQ.6) THEN 
WK33(3,3) = 1.D0 

END IF 


CALL MULTP2( BK ,WK33,16,3,3,WK163) 
CALL MULTP2(WK163,WK316,16,3,16,WK1616) 


C xk MULTIPLY WIGHT FC. AND SUM UP FOR ELEMENT STIFF. MATRIX 


DO 550 J=1, 16 





615 





DO 550 I=1, 16 
EK(I,J,KK) = EK(I,J,KK) + WT4(NI)*WT4(NJ)*WK1616(1,J) 
CONTINUE 
CONTINUE 


CONTINUE 


* ASSEMBLE FOR ELEMENT STIFFNESS MATRIX 


DO 620 I[=1, 6 
DO 615 J=1, 16 
DO 615 K=1, 16 
EK(K,J,1) = EK(K,J,1)*AAA*BBB 
CONTINUE 
DO 620 J=1, 16 
DO 620 K=1, 16 
EKT(K,J) = EXT(K,J) + DV(I)*EK(K,J,1) 
CONTINUE 


EI = EBEAM*BMINER 
ARSTF(1,1) = EI*12.D0/AAA**3 


ARSTF(J),2) = -~EI*6.D0/AAA**2 
ARSTF (1,3) = ~ARSTF(1,1) 
ARSTF(1,4) = ARSTF(1,2) 
ARSTF(2,2) = EI*4.D0/AAA 
ARSTF(2,3) = ~ARSTF(1,2) 
ARSTF(2,4) = EI*2.D0/AAA 
ARSTF(3,3) = ARSTF(1,1) 
ARSTF(3,4) = -ARSTF(1,2) 
ARSTF(4,4) = ARSTF(2,2) 
ARSTF(2,1) = ARSTF(1,2) 
ARSTF(3,1) = ARSTF(1,3) 
ARSTF(3,2) = ARSTF(2,3) 
ARSTF(4,1) = ARSTF(1,4) 
ARSTF(4,2) = ARSTF(2,4) 
ARSTF(4,3) = ARSTF(3,4) 


FACTOR = EBEAM*ABEAM/(EPLT*APLT + EBEAM*ABEAM ) 
ARGEM(1,1) = FACTOR*1.2D0/AAA 


ARGEM(1,2) = -FACTOR*0.1D0 
ARGEM(1,3) = -ARGEM(1,1) 
ARGEM(1,4) = ARGEM(1,2) 
ARGEM(2,2) = FACTOR*2.D0/15.DO*AAA 
ARGEM(2,3) = -ARGEM(1,2) 
ARGEM(2,4) = -FACTOR*AAA/30 .DO 
ARGEM(3,3) = ARGEM(1,1) 
ARGEM(3,4) = -ARGEM(1,2) 
ARGEM(4,4) = ARGEM(2,2) 
ARGEM(2,1) = ARGEM(1,2) 
ARGEM(3,1) = ARGEM(1,3) 
ARGEM(3,2) = ARGEM(2,3) 
ARGEM(4,1) = ARGEM(1,4) 
ARGEM(4,2) = ARGEM(2,4) 
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ARGEM(4,3) = ARGEM(3,4) 


RETURN 
END 


SUBROUTINE ASSEMB 


THIS SUBROUTINE ASSEMBLES THE ELEMENT STIFFNESS, ELEMENT GEOM. 
STIFFNESS MATRIX TO OBTAIN THE GLOBAL MATRIX 


IMPLICIT REAL*8(A-H,0-Z) 

COMMON/ELEMNT/ EK(16,16,6), EKT(16,16), EG(16,16) 

COMMON/ELEMNB/ ARSTF(4,4), ARGEM(4,4) 

COMMON/PROPT1/ E1(15), E2(15), G12(15), PO12(15), ANGL(15), 2(16) 
* » NOD(150,4), NX, NY, NSDF, LSDF(150) 

COMMON/PROPT3/ EBEAM, ABEAM, BMINER, EPLT, APLT, THICK, NODB(18,2) 
COMMON/TOTAL / GK(150,150) ,GG( 150,150) 

COMMON/BAND/ BANDK(45,150), BANDKR(45,150), BANDG(45,150) 
COMMON/EIG/NBAND, NREQ, NEQ 

DIMENSION NK(4) 


hoe 


CLEAR ARRAY 


DO 50 I=1, NEQ 
DO 50 J=1, NEQ 
GK(J,I) = 0.D0 
GG(J,I) = 0.D0 
CONTINUE 


NEL = NX*NY 
DO 100 N=1, NEL 


DO 100 I=1, 4 
NR=(NOD(N,I)-1)*4 
DO 100 II=1, 4 
NR=NR+1 
L=(I-1)*4+II 
DO 200 J=1, 4 
NCL=(NOD(N,J)-1)*4 
DO 200 JJ=1, 4 
M=(J-1)*4+JJ 
NC=NCL+JJ 
GK(NR ,NC)=GK(NR ,NC)+EKT(L,M) 
GG(NR ,NC)=GG(NR ,NC)+EG(L,M) 


200 CONTINUE 
100 CONTINUE 


C *** IMPOSITION OF THE BEAM ELEMENT MATRIX *** 
C 


DO 250 I=1, NX 
NOD1 = 4*(NODB(I,1)-1) 





NOD2 = 4*(NODB(I,2)-1) 


NK(1) = NOD1 + 1 
NK(2) = NOD1 + 3 
NK(3) = NOD2 + 1 
NK(4) = NOD2 + 3 


DO 240 J=1, 4 
DO 230 K=1, 4 
GK(NK(K) ,NK(J)) 


GG(NK(K) ,NK(J)) 
230 CONTINUE 
240 CONTINUE 
250 CONTINUE 
Cc 


C *** IMPOSITION OF THE BOUNDARY CONDITIONS *** 
Cc ‘ 
DO 310 NS=1,NSDF 
K=LSDF(NS)-NS+1 
DO 320 I=1,NEQ 
DO 320 J=K,NEQ 
GK(1I,J)=GK(I,J+1) 
GG(1I,J)=GG(I,J+1) 
320 CONTINUE 


DO 330 I=K,NEQ 
DO 330 J=1,NEQ 
GK(I,J)=GK(I+1,J) 
GG(I,J)=GG(I+1,J) 
330 CONTINUE 
310 CONTINUE 
C 
C **#* BAND STORAGE *** 
C 
NHBD = 4*(NX + 3) 
IF(NHBD.GE.NREQ) NHBD = NREQ 
NBAND = NHBD - 1 








C 
DO 20 J = 1, NREQ 
Il = MAX0(1, J-NBAND) 
po 10 I = 11, J 
K = I-J+NBAND+1 
BANDK(K,J) = GK(I,J) 
BANDKR(K,J) = GK(I,J) 
BANDG(K,J) = GG(I,J) 
10 CONTINUE 
20 CONTINUE 
C 
RETURN 
END 
Cc 
Cc 
C 
SUBROUTINE EIGEN 
C 
C 
Cc 


IMPLICIT REAL*8(A-H,0-Z) 
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GK(NK(K) ,NK(J)) + ARSTF(K,J) 
GG(NK(K) ,NK(J)) + ARGEM(K,J) 
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COMMON ENCE EGY E1(15), E2(15), G12(15), PO12(15), ANGL(15), 2(16) 
NOD(150,4), NX, NY, NSDF, LSDF(150) 

" counON/PROPTS) EBEAM, ABEAM, BMINER, EPLT, APLT, THICK, NODB(18,2) 

COMMON /PROPT4/ HT(15), AAA, BBB, NLYR, NLYR1, NLYR2, NLYR3 

COMMON/TOTAL / GK( 150,150) ,GG( 150, 150) 

COMMON/BAND/ BANDK(45,150), BANDKR(45,150), BANDG(45,150) 

COMMON/EIG/NBAND, NREQ, NEQ 

COMMON/EIG2/ EGVC(150,5), EGVCO(150,5), EGVL(5,4), EGVLO(S,4), 


* DTHK(50,5), DGD(5), EGVCB(150,5), IDMOD 


COMMON /NORMAL/ VNORMAL 


DIMENSION WORK(22500), IND(3) 
EQUIVALENCE (GG(1,1), WORK(1)) 
EXTERNAL OP, IOVECT 


FACTOR THE BANDED STIFFNESS MATRIX *** 
(LINPACK SUBROUTINE "“DPBFA") 


LDA = 45 
CALL DPBFA(BANDK, LDA, NREQ, NBAND, INFO) 
IFC INFO.NE.0) STOP 


INPUT FOR SUB-SPACE ITERATION *** 


NVAL = 3 
NFIG = 10 
NPERM = 0 
NMVAL = 3 
NMVEC = 150 


NBLOCK = NREQ/6 
IF(NBLOCK.GE.8) NBLOCK = 
MAXOP = NREQ 
MAXJ = 150 

NREQ2 = 2*NREQ 


CLEAR THE WORK SPACE 


NV = IABS(NVAL) 

NWORK1 = 2*NREQ*NBLOCK + MAXJ*(NBLOCK+NV+2) + 2*NBLOCK**2 +3*NV 
NWORK2 = N*NBLOCK 

NWORK3 = MAXJ*(2*NBLOCK+3) + 2*NV + 6 + (2*NBLOCK+2)*(NBLOCK+1) 
IF(NWORK2.GT.NWORK3) NWORK3 = NWORK2 

NWORK = NWORK1 + NWORK3 


INITIAL VALUE OF THE EIGENVECTORS 


NPR7 = NV 
IF(NBLOCK.LE.NV) NPR7 = NBLOCK 
DO 608 I=1, NPR7 
DO 607 J=1, NREQ 
WORK((I-1)*NREQ+J) = EGVCB(J,I) 
CONTINUE 
CONTINUE 


DO 1 I=1, NREQ 
WORK(I) = 0.D-05 
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1 CONTINUE 
C 
CALL DNLASO(OP, IOVECT, NREQ, NVAL, NFIG, NPERM, 
NMVAL, EGVL, NMVEC, EGVC, NBLOCK, MAXOP, MAXJ, 
WORK, IND, IERR) 
C 
PRINT 1400, IERR 
1400 FCRMAT( ' ERROR PARAMETER =',15) 
PRINT 1420, IND(1) 
1420 FORMAT(' NO. OF CALLS OF OP =',I5) 
PRINT 1430, NPERM 
1430 FORMAT(' NO. OF KNOWN EIGEN VALUES =',I5) 
IFC IERR.GT.0) STOP 
Cc 
WRITE(6,662) 
662 FORMAT(//,15X,'EIGEN SOLUTION FOR THE PLATE', 
/ 15X, 35('=") ) 
Cc 
C *** POSTPROCESS OF THE EIGENVALUES *** 
C 
DENO = EPLT*APLT + EBEAM*ABEAM 
DO 760 I=1, NV 
EGVL(I,1) = 1.D0/EGVL(I,1) 
SIG1 = EGVL(1I,1)*EPLT/DENO 
SIG2 = EGVL(I,1)*EJEAM/DENO 
Pl = SIG1*APLT 
P2 = SIG2*ABEAM 
VNX = SIG1*THICK 
Cc 
DO 250 I=1, NV 
Cress SAVE EGVCB 
DO 245 K=1, NREQ 
EGVCB(K,I) = EGVC(K,I) 
245 CONTINUE 
Cc 
C *** OBTAIN X USING THE SUBROUTINE "DPBSL2" *** 
Cc 
CALL DPBSL2(BANDK ,45 ,NRcQ,NBAND,®GVC(1.15))5 
Cc 
EGVL(1I,1) = EGVL(I,1)*VNORMAL 
WRITE(6,666) I, EGVL(I.1) 
250 CONTINUE 
RETURN 
END 
C 
c 
C 
Cc 
SUBROUTINE OP(N,M,P,Q) 
Cc 
C 


IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION P(N,M), Q(N,M) 

DIMENSION PBAR( 150% 

COMMON/BAND/ BANDK( 45,150), BANDAR(45 i50), BANDG(45,150) 
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COMMON/EIG/NBAND, NREQ, NEQ 


DO 999 L=1, M 
DO 105 J=1, N 
PBAR(J) = P(J,L) 
CONTINUE 


OBTAIN PBAR(N,M) USING THE SUBROUTINE "DPBSL2" *** 
CALL DPBSL2( BANDK ,45 ,NREQ,NBAND, PBAR(1) ) 

BANDED GEOMETRIC STIFFNESS MATRIX TIMES PBAR *** 

NBAND1 = NBAND + 1 


LC = 45 
CALL BNDMUL(BANDG, PBAR ,NBAND1,NREQ,Q(1,L),LC) ’ 


: OBTAIN Q(N,M) USING THE SUBROUTINE "DPBSL1" ** 


CALL DPBSL1(BANDK ,45 ,NREQ ,NBAND ,Q(1,L)) 
CONTINUE 
RETURN 
END 


Pi bee 


SUBROUTINE IOVECT(N,M,Q,J,K) 





IMPLICIT REAL*8(A~H,0-Z) 
COMMON/TOTAL / GK(150,150),GG(150,150) 
DIMENSION Q(N,M) 
IF(K.EQ.1) GO TO 30 
DO 20 L=1, M 

Li = J-M+tL 

DO 10 I=1, N 

GK(I,L1) = Q(I,L) 

CONTINUE 
CONTINUE 
RETURN 


DO 50 L=1, M 
Li = J-M+L 
DO 40 I=1, N 
Q(I,L) = GK(I,L1) 
CONTINUE 
CONTINUE 
RETURN 
END 


C 
OR se ee 
C 


Cc 
C 


SUBROUTINE BNDMUL( BAND ,WK ,NBAND1,NREQ,WK2,LC) 


LL 


IMPLICIT REAL*8(A-H,O-Z) 





4150 


4160 


4200 


agaaAaaAgaARAaANn aagaaa 


210 


220 
230 


DOUBLE PRECIGION BAND(LC,*), WK(*), WK2(*) 
NBAND = NBAND1 - 1 
DO 4200 I=1, NREQ 
SUM = 0.D0 
DO 4150 J=1, NBAND1 
IPR = I - NBAND1 + J 
IF(IPR.LE.0) GO TO 4150 
SUM = SUM + BAND(J,1I)*WKCIPR) 
CONTINUE 
DO 4160 K=1, NBAND 


IPR 


=I+kK 


IF(IPR.GT.NREQ) GO TO 4160 

SUM = SUM + BAND(NBAND1-K,IPR)*WK( IPR) 
CONTINUE 
WK2(1) = SUM 


CONTINUE 


RETURN 
END 


SUBROUTINE GRMESH(NNM, NEM) 


THIS SOUBROUTINE GENERATES THE BOOLEAN CONNECTIVITY MATRIX 


IMPLICIT REAL*8(A-H,0-Z2) 


COMMON/PROPT1/ E1(15), £2(€15), G12(15), PO12(15), ANGL(15), 2(16) 


* 


NX1 
NY1 


“wou 


NXX1=NXX+1 
NYY1=NYY+1 


NEM 


, NOD(150,4), NX, NY, NSDF, LSDF(150) 


NX +1 
NY +1 


= NX*NY 


NNM = NX1*NY1 


NOD(1,1) 
NOD(1,2) 
NOD(1,3) 
NOD(1,4) 


1 
2 
NX1 + 2 
NX1 + 1 


IF(NY.EQ.1, GO TO 230 


M=1 


DO 220 N=2, NY 
L = (N-1)*NX + 1 
DO 210 I=1, 4 
NOD(L,I) = NOD(M,I) + NX1 
CONTINUE 


M=L 


CONTINUE 


IF(NX.EQ.1) GO TO 270 
DO 260 NI=2, NX 
DO 240 I=1, 4 


7] 





ak TF] 


240 


250 


260 
270 


NOD(NI,I) = NOD(NI-1,1) + 1 


CONTINUE 
M=NI 
DO 260 NJ=2, NY 
L = (NJ-1)*NX + NI 
DO 250 J=1, 4 
NOD(L,J) = NOD(M,J) + NX1 
CONTINUE 
M=L 
CONTINUE 
RETURN 
END 


ieee 


aaAaAaaANn AQ 


aQaaa 


aaa 


aaa 





APPENDIX B. OPTIMIZATION PROBLEM CODE LISTING USING 
IMSL 


PROGRAM MAINIM 


VARIABLE DECLARATIONS 


IMPLICIT REAL*8(A~-H,0-Z) 

DOUBLEPRECISION XLB(18), XUB(18), XGUESS(18), T(18) 
PARAMETER (IBTYPE=0, IPRINT=3, M=5, MAXITN=100, ME=1) 
EXTERNAL FCN, GRAD, INPUTB, LIMITS, GUESS, DNCONG 


COMMON BLOCKS 


COMMON/ELEMNT/ EK(16,16,6), EKT(16,16), EG(16,16) 
COMMON/ELEMNB/ ARSTF(4,4), ARGEM(4,4) 
COMMON/PROPTO/ THETA, ALEN, BLEN 
COMMON/PROPT1/ E1(15), E2(15), G12(15), PO12(15), ANGL(15), 2(16), = 
2 NOD(150,4), NX, NY, NSDF, LSDF(150) z 
COMMON/PROPT2/ DM(32,3) 
COMMON/PROPT3/ EBEAM, ABEAM, BMINER, EPLT, APLT, THICK, NODB(18,2) 
. COMMON/PROPT4/ HT(15), AAA, BBB, NLYR, NLYR1, NLYR2, NLYR3 
COMMON/TOTAL/ GK(150,150), GG( 150,150) 
COMMON/BAND/ BANDK(45,150), BANDKR(45,150), BANDG(45,150) 
COMMON/EIG/ NBAND, NREQ, NEQ 
COMMON/EIG2/ EGVC(150,5), EGVCO(150,5), EGVL(5,4), EGVLO(5,4), 
2 DTHK(50,5), DGD(5), EGVCB(150,5), IDMOD 
COMMON/FSAVE/ VVV2(100) 
COMMON/BEAM/ £E1B, E2B, G12B, PO12B, BBM, HBM 
COMMON/BEAM2/ DBTHK(50), PBEAM 
COMMON/BIMO/ VECRAT(S,5), V99, INDVC(5) 
COMMON/LEVY/ SD11, SD12, SD22, SD66, ALPM, ALPMS, PEM, OMEGA 
COMMON/NORMAL/ VNORMAL 


VNORMAL IS SCALING FACTOR FOR T(NLYR3) TO INCREASE ACCURACY 
VNORMAL = 1.0D 03 


CALL INPUTB 

N = NLYR + 3 

CALL LIMITS (XLB, XUB) 

CALL GUESS (XGUESS) 

CALL DNCONG (FCN, GRAD, M, ME, N, XGUESS, IBTYPE, XLB, XUB, 
2 IPRINT, MAXITN, T, PMAX) 

END 
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Cc 
C 
C 
C 


C 





SUBROUTINE FCN (M, ME, N, T, ACTIVE, P, G) 





IMPLICIT REAL*8(A-H,0-Z) 
LOGICAL ACTIVE(*) 
DOUBLEPRECISION G(*), T(*) 
EXTERNAL INPUTC, EIGANL4, LEVYS 


COMMON BLOCKS *xxterteictkaitketkteee eke 


COMMON /ELEMNT/ 
COMMON /ELEMNB/ 
COMMON/PROPTO/ 
COMMON /PROPT1/ 


COMMON/PROPT2/ 
COMMON / PROPT3/ 
COMMON/PROPT4/ 
COMMON /TOTAL/ 
COMMON /BAND/ 
COMMON/EIG/ 
COMMON/EIG2/ 


2 


COMMON /BEAM/ 
COMMON /BEAM2/ 


EK(16,16,6), EKT(16,16), EG(16,16) 
ARSTF(4,4), ARGEM(4,4) 
THETA, ALEN, BLEN 


E£1(15), E2(15), G12(15), PO12(15), ANGL(15), 2(16), : 
NOD(150,4), NX, NY, NSDF, LSDF(150) 

DM(3,3) 

E2EAM, ABEAM, BMINER, EPLT, APLT, THICK, NODB(18,2) _ 


i. 15), AAA, BBB, NLYR, NLYR1, NLYR2, NLYR3 

Gx( 150,150), GG(150,150) 

BANDK(45,150), BANDKR(45,150), BANDG(45,150) 
NBAND, NREQ, NEQ 

EGVC(150,5), EGVCO(150,5), EGVL(5,4), EGVLO(5,4), 
DTHK(50,5), DGD(5), EGVCB(150,5), IDMOD 

E1B, E2B, G12B, PO12B, BBM, HBM 

DBTHK(50), PBEAM 


 fpre 


C UPDATE VALUES OF HT(*), BBM, AND HBM 


C 


aaaa 


aan aaa 


DO 100 I=1,NLYR 


BBM 
HBM 


CALCUL 


HT(I) = TCT) 
100 CONTINUE 


TCNLYR1) 
TCNLYR2) 


TE VALUES OF Pl, P2, P3, AND PBEAM FOR CURRENT VALUES OF 
HT(*), BBM, AND HBA 


INPUTC 
EIGANL4 
LEVYS 
EGVL(1,1) 
EGVL(2,1) 
EGVL(3,1) 


FUNCTION DEFINITION 


TCNLYR3) 


CALCULATE CONSTRAINT VALUES FOR CURRENT VALUES OF HT(*), BBM, & HBM 


IF (ACTIVE(1)) G(1) 
IF (ACTIVE(2)) G(2) 
IF (ACTIVE(3)) G(3) 
IF (ACTIVE(4)) G(4) 





APLT + ABEAM - THETA 
Pl - TC(NLYR3) 
.999*P2 - T(NLYR3) 
.998*P3 - T(NLYR3) 
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aQaaqang aa aqaaa aQ 


agaaqa anaaq a ann 


aaa 


Cc 


IF (ACTIVE(S)) G(5) = PBEAM - T(NLYR3) 


RETURN 


SUBROUTINE GUESS (XGUESS) 


IMPLICIT REAL*8(A-H,0-Z) 
DOUBLEPRECISION XGUESS( 18) 


COMMON/PROPT4/ HT(15), AAA, BBB, NLYR, NLYR1, NLYR2, NLYR3 


COMMON/NORMAL/ VNORMAL 


READ IN XGUESS VALUES 


DO 200 I=1,NLYR3 
READ(8,*) XGUESS(I) 


200 CONTINUE 


SCALE T(NLYR3) GUESS FOR ADDITIONAL ACCURACY 


XGUESS (NLYR3 )=XGUESS (NLYR3 )*VNORMAL 
RETURN 


SUBROUTINE LIMITS (XLB, XUB) 





IMPLICIT REAL*8(A-H,0-Z) 

COMMON/PROPTO/ THETA, ALEN, BLEN 

COMMON/PROPT4/ HT(15), AAA, BBB, NLYR, NLYR1, NLYR2, 
DOUBLEPRECISION XLB(18), XUB(18) 


SET LOWER LIMITS FOR DESIGN VARIABLES, T(*) 


100 


DO 100 I=1,NLYR 
XLB(I) = 1.0D-03 
CONTINUE 
XLB(NLYR1) 
XLB(NLYR2 ) 
XLB(NLYR3) 


1.0D-03 
1.0D-03 
1.0D-06 


C SET UPPER LIMITS FOR DESIGN VARIABLES, HT(*) 


C 


200 


DO 200 I=1,NLYR 
XUB(I) = 2.0D-01 
CONTINUE 
XUB(NLYR1) 
XUB(NLYR2) 
XUB(NLYR3) 


nou 
hm KO DO 
ooo 
oO 

' 

oO 

i) 


RETURN 
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NLYR3 





APPENDIX C. OPTIMIZATION PROBLEM CODE LISTING USING 


DOT 
Cc 
C 
PROGRAM MNDOT 
Cc 
Cc 
C VARIABLE DECLARATIONS 
Cc 
IMPLICIT REAL*8(A-H,0-Z) ‘ 
DOUBLEPRECISION XLB(18), XUB(18), XGUESS(18), T(18), DG(6,18), 
2 G(6), AAC18,6), WK(2000), RPRM(20) 
DIMENSION IPRM(20), IWK(1000) ; 
EXTERNAL FCNDOT, INPUTB, LIMITS, GUESS, DOT 
Cc 
C COMMON BLOCKS 
c 
COMMON/ELEMNT/ EK(16,16,6), EKT( 16,16), EG(16,16) 
COMMON/ELEMNB/ ARSTF(4,4), ARGEM(4,4) 
COMMON/PROPTO/ THETA, ALEN, BLEN = 
COMMON/PROPT1/ E1(15), E2(15), G12(15), PO12(15), ANGL(15), Z(16), & 
2 NOD(150,4), NX, NY, NSDF, LSDF(150) @ 
COMMON/PROPT2/ DM(3,3) 
COMMON/PROPT3/ EBEAM, ABEAM, BMINER, EPLT, APLT, THICK, NODB(18,2) 
COMMON/PROPT4/ HT(15), AAA, BBB, NLYR, NLYR1, NLYR2, NLYR3 
COMMON/TOTAL/ GK(150,150), GG( 150,150) 
COMMON/BAND/ BANDK(45,150), BANDKR(45,150), BANDG(45,150) 
COMMON/EIG/ NBAND, NREQ, NEQ 
COMMON/EIG2/ EGVC(150,5), EGVCO(150,5), EGVL(5,4), EGVLO(5,4), 
2 DTHK(50,5), DGD(5), EGVCB(150,5), IDMOD 
COMMON/FSAVE/ VVV2(100) 
COMMON/BEAM/ £E1B, E2B, G12B, PO12B, BBM, HBM 
COMMON/BEAM2/ DBTHK(50), PBEAM 
COMMON/BIMO/ VECRAT(S,5), V99, INDVC(5) 
COMMON/LEVY/ SD11, SD12, SD22, SD66, ALPM, ALPMS, PEM, OMEGA 
COMMON/NORMAL/ VNORMAL, SCALTH 
Cc 
C VNORMAL IS SCALING FACTOR FOR T(NLYR3) TO INCREASE ACCURACY ‘ 
Cc 
VNORMAL = 1.0D+02 
* 
Cc 
C DEFINE WORK SPACE REQUIREMENTS 
Cc 
NRWK = 2000 
NRIWK = 1000 
C 
C ZERO RPRM AND IPRM 
Cc 


DO 50 I[=1,20 
RPRM(I) = 0.0D0 





IPRM(I) = 0 
50 CONTINUE 


C 
C READ IN PROBLEM DATA AND VALUES OF T(*), XLB(*), AND XUB(*) 
C 
CALL INPUTB 
CALL LIMITS (XLB, XUB) 
CALL GUESS (XGUESS) 
DO 100 I=1,NLYR3 
TCI) = XGUESS(I) 
100 CONTINUE 





Cc 
C SPECIFY THAT GRADIENTS ARE TO BE PROVIDED 
C SET PARAMETER IPRM(1) = 1 
Cc 
IPRM(1) = 1 
C 
C DEFINE METHOD, NDV, NCON, IPRINT, MINMAX 
C METHOD = 1 FOR MFD 
C METHOD = 2 FOR SLP 
C 
METHOD = 
NDV = NLYR3 
NCON = 5 
IPRINT = 5 - 
MINMAX = 0 Fe 
sl 
C READY TO OPTIMIZE 
Cc 
INFO = 0 
200 CALL DOT (INFO,METHOD,IPRINT,NDV,NCON,T,XLB,XUB,P,MINMAX,G, 
2 RPRM, IPRM,WK ,NRWK, IWK,NRIWK) 
a 
C PROVIDE GRADIENTS IF DOT REQUESTS THEM 
Cc 
IFC INFO.EQ.2) GO TO 300 
Cc 
C EXIT IF CONVERGENCE IS OBTAINED 
Cc 
IF(INFO.EQ.0) GO TO 1000 
Cc 
C EVALUATE OBJECTIVE FUNCTION AND CONSTRAINTS 
Cc 
CALL FCNDOT (T, P, G) 
GO TO 200 
Cc 
C GRADIENT EVALUATION 
Cc 
300 CONTINUE 
Cc 


C CALCULATE FUNCTION GRADIENT WRT DESIGN VARIABLES AT POINT HT(*) 
C 


DO 400 I=1,NI-YR2 
WK(I) = 0.0D0 
400 CONTINUE 





aaa 


C 
C 
Cc 


WK(NLYR3) = -1.0D0 


IF NUMBER OF ACTIVE CONSTRAINTS EQUAL ZERO RETURN TO DOT 


NGT = IPRM(20) 


IF(NGT.EQ.0) GO TO 200 


CALCULATE CONSTRAINT GRADIENTS WRT DESIGN VARIABLES AT POINT HT(*) 


DO 510 I=1,NLYR 


DG(1,I) = -DTHK(I,1) 


STORE ACTIVE GRADIENT INFORMATION IN ARRAY WK(*) 


510 CONTINUE 

DG(1,NLYR3) = 1.0D0 
DO 520 I=1,NLYR2 
DG(2,1) = -DTHK(I,2) * 0.999D0 

520 CONTINUE 

DG(2,NLYR3) = 1.0D0 
DO 530 I=1,NLYR2 
DG(3,I) = -DTHK(I,3) * 0.998D0 

530 CONTINUE 

DG(3,NLYR3) = 1.0D0 
DO 540 I=1,NLYR2 
DG(4,I1) = -DBTHK(I) 

540 CONTINUE 

DG(4,NLYR3) = 1.0D0 
DO 550 I=1,NLYR 
DG(5,I) = 2. 

550 CONTINUE 
DG(5,NLYR1) = 
DG(S,NLYR2) = 
DG(5,NLYR3) = 0.0D0 

DO 600 K=1,NGT 
J = IWK(K) 
AA(1,K) = DG(J,1) 
AA(2,K) = DG(J,2) 
AA(3,K) = DG(J,3) 
AA(4,K) = DG(J,4) 
AA(S,K) = DG(J,5) 
600 CONTINUE 
N1 = NDV 
DO 800 K=1,NGT 
DO 700 I=1,NDV 
WK(I+N1) = AA(I,K) 
700 CONTINUE 


Nl = Nl + NDV 


800 CONTINUE 
GO TO 200 
1000 CONTINUE 


2 


ODO * BLEN/THETA 


2.0D0 * T(NLYR2)/THETA 
2.0D0 * T(NLYR1)/THETA 
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Pipes 


ana aa 


aan 


Cc 


END 








SUBROUTINE FCNDOT (T, P, G) 


IMPLICIT REAL*8(A-H,0-Z) 
DOUBLEPRECISION G(*), T(*) 
EXTERNAL INPUTC, EIGANL4, LEVYS 


COMMON BLOCKS ***xtatrkkitkkkkekiiikiii 


COMMON /ELEMNT/ 


COMMON /ELEMNB/ 
COMMON /PROPTO/ 
COMMON /PROPT1/ 


2 


COMMON / PROPT2/ 
COMMON /PROPT3/ 
COMMON/PROPT4/ 
COMMON/TOTAL/ 
COMMON /BAND/ 
COMMON /EIG/ 
COMMON/EIG2/ 


2 


COMMON / BEAM / 


COMMON/BEAM2/ DBTHK(50), PBEAM 


EK(16,16,6), EKT(16,16), EG(16,16) 

ARSTF(4,4), ARGEM(4,4) 

THETA, ALEN, BLEN 

E1(15), £2(€15), G12(15), PO12(15), ANGL(15), 2(16), 
NOD(150,4), NX, NY, NSDF, LSDF(150) 

DM(3,3) 

EBEAM, ABEAM, BMINER, EPLT, APLT, THICK, NODB(18,2) 
HT(15), AAA, BBB, NLYR, NLYR1, NLYR2, NLYR3 
GK(150,150), GG(150,150) 

BANDK(45,150), BANDKR(45,150), BANDG(45,150) 

NBAND, NREQ, NEQ 

EGVC(150,5), EGVCO(150,5), EGVL(5,4), EGVLO(5,4), 
DTHK(50,5), DGD(5), EGVCB(150,5), IDMOD 

E1B, E2B, G12B, PO12B, BBM, HBM 


a de 


COMMON/NORMAL/ VNORMAL, SCALTH 


C UPDATE VALUES OF HT(*), BBM, AND HBM 


Cc 


aaAaaAn 


aaa aaa 


100 


DO 100 I=1,NLYR 
HT(I) = TCI) 
CONTINUE 


BBM 
HBM 


T(NLYR1) 
TC(NLYR2) 


CALCULATE VALUES OF P1, P2, P3, AND PBEAM FOR CURRENT VALUES OF 
HT(*), BBM, AND HBM 


CALL 
CALL 
CALL 
Pl = 
P2 
P3 


INPUTC 
EIGANL4 
LEVYS 
EGVL(1,1) 
EGVL(2,1) 
EGVL(3,1) 


CALCULATE FUNCTION VALUE FOR CURRENT VALUES OF HT(*), BBM, & HBM 


P = -T(NLYR3) 


CALCULATE 


G(1) 
G(2) 
G(3) 


CONSTRAINT VALUES FOR CURRENT VALUES OF HT(*), BBM, & HBM 


(T(NLYR3)-P1) 
(T(NLYR3)-0.999D0*P2) 
(TCNLYR3)-0.998D0*P3) 
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G(4) 
G(5) 


(TC(NLYR3)-PBEAM ) 
(APLT + ABEAM)/THETA - 1.0D0 


RETURN 


aaa Q 


SUBROUTINE GUESS (XGUESS) 


aa 


IMPLICIT REAL*8(A-H,0-Z) 
DOUBLEPRECISION XGUESS( 18) 


COMMON/PROPT4/ HT(15), AAA, BBB, NLYR, NLYR1, NLYR2, NLYR3 


COMMON/NORMAL/ VNORMAL 
READ IN XGUESS VALUES 


qaqa 


DO 200 I=1,NLYR3 
READ(8,*) XGUESS(I) 
200 CONTINUE 
SCALE T(NLYR3) GUESS FOR ADDITIONAL ACCURACY 
XGUESS(NLYR3 )=XGUESS(NLYR3)*VNORMAL 


RETURN 





SUBROUTINE LIMITS (XLB, XUB) 





aa aaa Q aaa 


IMPLICIT REAL*8(A~-H,0-Z) 
COMMON/PROPTO/ THETA, ALEN, BLEN 


COMMON/PROPT4/ HT(15), AAA, BBB, NLYR, NLYR1, NLYR2, NLYR3 


DOUBLEPRECISION XLB(18), XUB(18) 
SET LOWER LIMITS FOR DESIGN VARIABLES, T(*) 


aqaa 


DO 100 I=1,NLYR 
XLB(I) = 1.0D-03 
100 CONTINUE 


XLBC(NLYR1) = 1.0D-03 
XLB(NLYR2) = 1.0D-03 
XLB(NLYR3) = 1.0D-06 


C 
C SET UPPER LIMITS FOR DESIGN VARIABLES, HT(*) 
C 
DO 200 I=1,NLYR 
XUB(I) = 2.0D-01 
200 CONTINUE 


XUB(NLYR1) = 2.0D-01 

XUB(NLYR2) = 2.0D-01 

XUBC(NLYR3) = 1.0D 06 
C 

RETURN 
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